diff --git a/README.md b/README.md index e39f9a5..00b523b 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,6 @@ # assignment-bank-2015 -For submission of class projects at GW, Fall 2015 +For submission of class project at GW, Fall 2015 + +This assignment examines numerical integration techniques using mechanical systems. In the notebook Integration_Techniques_for_Mechanical_Systems, I start by introducing a simple harmonic oscillator and three different types of integrators (Euler, Runge-Kutta, and Verlet). Next I introduce damping and look at how the system behavior changes for each method of integration. The Double_Pendulum_Problem notebook looks at the double pendulum problem using Runge-Kutta and Verlet integration. + +The code in this repository is under Creative Commons Attribution license CC-BY 4.0, code under MIT license © 2015 Randall Schur diff --git a/randy_schur/Double_Pendulum_Problem.ipynb b/randy_schur/Double_Pendulum_Problem.ipynb new file mode 100644 index 0000000..e9c8071 --- /dev/null +++ b/randy_schur/Double_Pendulum_Problem.ipynb @@ -0,0 +1,7181 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Dynamics of a Double Pendulum" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this notebook, we will study the dynamics of a double pendulum. Like before, we'll derive the equations of motion from the Lagrangian and the Hamiltonian. We can use the same methods of integration as before, but this time there is no analytical solution to compare to. Finally, we will look at the properties of the system in phase space using Poincare sections." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Background\n", + "\n", + "The double pendulum problem is a relatively simple system which can produce surprisingly complex motion. It is a chaotic system, meaning it is unpredictable and small changes in initial conditions lead to large changes in motion. There is no closed form solution to the motion of the two masses. This is a frequent example problem in the study of dynamic systems, nonlinear controls, and mechanics.\n", + "\n", + "In this notebook we will be treating the two rods as massless and we'll start out ignoring the effects of friction. For a treatment of this problem without these assumptions, see [4].\n", + "\n", + "See the image below for the definition of constants." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "![Image](figures/diagram.png)\n", + "#### Figure 1. Diagram of Double Pendulum System" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Equations of Motion \n", + "\n", + "As a reminder, the Lagrangian is\n", + "$$\\begin{equation*}\n", + "L = T- U\n", + "\\end{equation*}$$\n", + "where $T$ is the total Kinetic Energy in the system, and $U$ is the total Potential Energy. Potential energy is calculated as:\n", + "\n", + "$$\\begin{eqnarray*}\n", + "U &=& m_1g(-y_1) + m_2g(-y_2) \\\\\n", + " &=& -m_1gl_1cos\\theta_1 - m_2g(l_1cos\\theta_1 + l_2cos\\theta_2)\n", + "\\end{eqnarray*}$$\n", + "\n", + "Kinetic Energy is calculated as:\n", + "\n", + "$$\\begin{eqnarray*}\n", + "T &=& \\frac{1}{2}mv_1^2+\\frac{1}{2}mv_2^2 \\\\\n", + "&=& \\frac{1}{2}m_1(l_1\\dot\\theta_1)^2 + \\frac{1}{2}m_2(l_1\\dot\\theta_1 + l_2\\dot\\theta_2)^2 \\\\\n", + "&=& \\frac{1}{2}(m_1+m_2)l_1^2\\dot\\theta_1^2+\\frac{1}{2}m_2\\left(2l_1l_2\\dot\\theta_1\\dot\\theta_2cos(\\theta_2-\\theta_1)+l_2^2\\dot\\theta_2^2\\right)\n", + "\\end{eqnarray*}$$\n", + "\n", + "\n", + "So the Lagrangian quantity becomes:\n", + "$$\\begin{equation*}\n", + "L = \\frac{1}{2}(m_1+m_2)l_1^2\\dot\\theta_1^2 + \\frac{1}{2}m_2l_2^2\\dot{\\theta_2}^2+m_2l_1l_2\\dot\\theta_1\\dot\\theta_2cos(\\theta_2-\\theta_1)+(m_1+m_2)l_1gcos\\theta_1 + m_2l_2gcos\\theta_2\n", + "\\end{equation*}$$\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "The equations of motion are then found using the Lagrange Equation:\n", + "\n", + "$$\\begin{equation*}\n", + "\\frac{d}{dt} \\left( \\frac{\\partial L}{\\partial \\dot{x_i}}\\right) - \\frac{\\partial L}{\\partial x_i} = 0\n", + "\\end{equation*}$$\n", + "\n", + "\n", + "There are two generalized coordinates ($\\theta_1$ and $\\theta_2$). So the equation of motion for $\\theta_1$ is calculated by the following steps:\n", + "\n", + "$$\\begin{eqnarray*}\n", + "\\frac{\\partial L}{\\partial \\dot{\\theta_1}} &=& (m_1+m_2)l_1^2\\dot\\theta_1+m_2l_1l_2\\dot\\theta_2cos(\\theta_2-\\theta_1) \\\\\n", + "\\frac{d}{dt} \\left( \\frac{\\partial L}{\\partial \\dot{\\theta_1}}\\right) &=& (m_1+m_2)l_1^2\\ddot\\theta_1+m_2l_1l_2\\ddot\\theta_2cos(\\theta_2-\\theta_1)-m_2l_1l_2\\dot\\theta_2^2sin(\\theta_2-\\theta_1)+m_2l_1l_2\\dot\\theta_1\\dot\\theta_2sin(\\theta_2-\\theta_1) \\\\\n", + "-\\frac{\\partial L}{\\partial \\theta_1} &=& -m_2l_1l_2\\dot\\theta_1\\dot\\theta_2sin(\\theta_2-\\theta_1)+(m_1+m_2)gl_1sin\\theta_1\n", + "\\end{eqnarray*}$$\n", + "\n", + "We can add together the second and third equations to find the first equation of motion. Some of these terms will cancel out, and we are left with a final equation of motion:\n", + "\n", + "$$\\begin{equation*}\n", + "0 = (m_1+m_2)l_1^2\\ddot\\theta_1+m_2l_1l_2\\ddot\\theta_2cos(\\theta_2-\\theta_1)-m_2l_1l_2\\dot\\theta_2^2sin(\\theta_2-\\theta_1)+(m_1+m_2)gl_1sin\\theta_1\n", + "\\end{equation*}$$\n", + "\n", + "We will also need the equation of motion for $\\theta_2$. This can be derived using the same process. It should work out to the following (check this by yourself!):\n", + "\n", + "$$\\begin{equation*}\n", + "0 = m_2l_2^2\\ddot\\theta_2+m_2l_1l_2\\ddot\\theta_1cos(\\theta_2-\\theta_1) + m_2l_2l_1\\dot\\theta_1^2sin(\\theta_2-\\theta_1)+ l_2m_2gsin\\theta_2\n", + "\\end{equation*}$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Problem setup\n", + "\n", + "Notice that these equations of motion are implicit equations for $\\ddot\\theta_1$ and $\\ddot\\theta_2$. Let's rewrite the equations of motion in a more readable form. We can use the standard form of this equation (which you'll be familiar with if you have studied any dynamics) by separating out mass terms (anything multiplied by $\\ddot{x}$), centripetal terms (anything multiplied by $\\dot{x}$), and gravitational terms (anything multiplied by $x$). We get the following equation of motion:\n", + "\n", + "$$\\begin{eqnarray*}\n", + "0 &=&\\textbf{M}\\vec{\\ddot{x}}(t) + \\textbf{V}\\dot{x} + \\textbf{G} \n", + "\\end{eqnarray*}$$\n", + "\n", + "The equation of motion can be solved for $\\ddot{x}$:\n", + "\n", + "$$ \\ddot{x} = \\textbf{M}^{-1}\\left[ \\textbf{V}\\dot x + \\textbf{G} \\right]$$\n", + "\n", + "where \n", + "\n", + "\\begin{eqnarray*}\n", + "M = \n", + "\\left[\\begin{array}{c}\n", + "(m_1+m_2)l_1 & m_2l_1l_2cos(\\theta_2 - \\theta_1)\\\\\n", + " m_2l_1l_2cos(\\theta_2 - \\theta_1) &m_2l_2^2\n", + "\\end{array}\\right]\n", + "\\end{eqnarray*}\n", + "\n", + "\n", + "\\begin{equation*}\n", + "V = \n", + "\\left[\\begin{array}{c}\n", + "0 & -m_2l_1l_2sin(\\theta_2 - \\theta_1)\\\\\n", + " m_2l_1l_2sin(\\theta_2 - \\theta_1) &0\n", + "\\end{array}\\right]\n", + "\\end{equation*}\n", + "\n", + "\n", + "\\begin{eqnarray*}\n", + "G = \n", + "\\left[\\begin{array}{c}\n", + "(m_1+m_2)l_1gsin\\theta_1\\\\\n", + " m_2l_2gsin\\theta_2\n", + "\\end{array}\\right]\n", + "\\end{eqnarray*}\n", + "\n", + "and\n", + "\n", + "$$\\begin{eqnarray*}\n", + "\\vec{x}(t) = \\begin{bmatrix} \\theta_1\\\\ \\theta_2\\\\ \\end{bmatrix}, \\vec{\\dot{x}(t)} = \\begin{bmatrix} \\dot \\theta_1 \\\\ \\dot \\theta_2 \\end{bmatrix}\n", + "\\end{eqnarray*}$$\n", + "\n", + "We can now treat this system the same way we handled the simple harmonic motion: break the two second order ODEs into four first-order ODEs, then use Euler or Runge-Kutta on the results." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Equations of Motion from Hamiltonian\n", + "\n", + "Here we derive the equations of motion from the Hamiltonian. If you want good practice, you can derive these on your own. If you want to see another derivation with nice animations, you can see [1]. We start with the same potential and kinetic energy as before, but we'll need to use postion ($\\theta$) and momentum ($I\\omega$) to give us:\n", + "\n", + "$$\\begin{eqnarray*}\n", + "U &=& m_1g(-y_1) + m_2g(-y_2) \\\\\n", + " &=& -m_1gl_1cos\\theta_1 - m_2g(l_1cos\\theta_1 + l_2cos\\theta_2) \\\\\n", + " &=& -m_1gl_1cos(q_1) - m_2g(l_1cos(q_1) + l_2cos(q_2))\n", + "\\end{eqnarray*}$$\n", + "\n", + "Kinetic Energy for rotational motion calculated as:\n", + "\n", + "$$\\begin{eqnarray*}\n", + "T &=& \\frac{1}{2}(m_1+m_2)l_1^2\\dot\\theta_1^2+\\frac{1}{2}m_2\\left(2l_1l_2\\dot\\theta_1\\dot\\theta_2cos(\\theta_2-\\theta_1)+l_2^2\\dot\\theta_2^2\\right)\\\\\n", + "&=& \\frac{1}{2}(I_1 + m_2l_1^2)\\omega_1^2+\\left[m_2l_1l_2\\omega_1\\omega_2cos(\\theta_2-\\theta_1)+\\frac{I_2\\omega_2^2}{2}\\right] \\\\\n", + "&=& \\frac{p_1^2}{2I_1}+\\frac{m_2p_1^2}{2m_1I_1} + \\left[ \\frac{p_2}{l_2} \\frac{p_1}{m_1l_1}cos(q_2-q_1) + \\frac{p_2^2}{2I_2} \\right]\n", + "\\end{eqnarray*}$$\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $I_i$ is the angular moment of inertia and $\\omega$ is the angular velocity. So our Hamiltonian quantity becomes:\n", + "$$\\begin{equation*}\n", + "H = \\frac{p_1^2}{2I_1}+\\frac{m_2p_1^2}{2m_1I_1} + \\left[ \\frac{p_2}{l_2} \\frac{p_1}{m_1l_1}cos(q_2-q_1) + \\frac{p_2^2}{2I_2} \\right] -m_1gl_1cos(q_1) - m_2g(l_1cos(q_1) + l_2cos(q_2))\n", + "\\end{equation*}$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Since we have two degrees of freedom, we will have four equations of motion. These will be the following:\n", + "$$\\begin{eqnarray*}\n", + "\\dot{q_1} &=& -\\frac{\\partial H}{\\partial p_1} = -\\frac{p_1}{I_1} - \\frac{m_2p_1}{m_1I_1} - \\frac{p_2}{m_1l_1l_2}cos(q_2-q_1) \\\\\n", + "\\dot{p_1} &=& \\frac{\\partial H}{\\partial q_1} = \\frac{p_1p_2}{m_1l_1l_2}cos(q_2-q_1) + m_1gl_1sin(q_1)+ m_2gl_2sin(q_1)\\\\\n", + "\\dot{q_2} &=& -\\frac{\\partial H}{\\partial p_2} = -\\frac{p_1}{m_1l_1l_2}cos(q_2-q_1) - \\frac{p_2}{I_2}\\\\\n", + "\\dot{p_2} &=& \\frac{\\partial H}{\\partial q_2} = -\\frac{p_1p_2}{m_1l_1l_2}sin(q_2-q_1) + m_2gl_2sin(q_2)\\\\\n", + "\\end{eqnarray*}$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Integration\n", + "\n", + "For Euler and Runge-Kutta, we can use the matrix equation above to integrate our system. The results will be VERY dependent on initial conditions; this is why we call it a chaotic system. For the symplectic integrator, we'll use the Hamiltonian equations of motion above. " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import numpy\n", + "from scipy.linalg import solve\n", + "from numpy.linalg import det\n", + "from math import pi, cos, sin, sqrt\n", + "\n", + "from matplotlib import pyplot\n", + "from matplotlib.pyplot import quiver\n", + "%matplotlib notebook\n", + "from matplotlib import rcParams, cm\n", + "rcParams['font.family'] = 'serif'\n", + "rcParams['font.size'] = 16" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def f_double_pendulum(u):\n", + " \"\"\"Returns RHS of double pendulum EOM\n", + " \n", + " Parameters:\n", + " u - initial state\n", + " \n", + " Returns:\n", + " RHS - RHS of harmonic oscillator eqn.\n", + " \n", + " \"\"\"\n", + " x1 = u[0]\n", + " x2 = u[1]\n", + " x3 = u[2]\n", + " x4 = u[3]\n", + " \n", + " M = numpy.array([[(m1+m2)*l1**2, m2*l1*l2*cos(x1-x2)],[m2*l1*l2*cos(x1-x2), m2*l2**2]])\n", + " V = numpy.array([[0, m2*l1*l2*sin(x1-x2)],[-m2*l1*l2*sin(x1-x2), 0]])\n", + " G = numpy.array([[(m1+m2)*g*l1*sin(x1)],[m2*l2*g*sin(x2)]])\n", + " qdd = numpy.linalg.inv(M).dot(V.dot(numpy.array([[x3],[x4]])) +G) #- 2*numpy.array([[x3],[x4]])) \n", + " RHS = numpy.array([x3, x4, qdd[0], qdd[1]])\n", + " return RHS" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "#Set up parameters:\n", + "g = -9.8 #[m/s^2]\n", + "m1 = 2\n", + "m2 = 1.5;\n", + "l1 = 1\n", + "l2 = 1.5\n", + "I1 = m1*l1**2\n", + "I2 = m2*l2**2\n", + "\n", + "T = 50; #[seconds]\n", + "dt = .01; #\n", + "N = int(T/dt)+1\n", + "t = numpy.linspace(0.0, T, N)\n", + "\n", + "#Initial Conditions\n", + "theta1_0 = pi/2 #[radians]\n", + "theta2_0 = pi\n", + "theta1_dot_0 = 0 #[rad/s]\n", + "theta2_dot_0 = 0 #\n", + "\n", + "q1_0 = theta1_0\n", + "p1_0 = m1*theta1_dot_0\n", + "q2_0 = theta2_0\n", + "p2_0 = m2*theta2_dot_0\n", + "\n", + "x_init_dp = numpy.array([theta1_0, theta2_0, theta1_dot_0, theta2_dot_0])\n", + "x_init_H_dp = numpy.array([q1_0, p1_0, q2_0, p2_0])" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def euler_DP(u, f, dt):\n", + " \"\"\" Euler's method for integrating a system of differential equations.\n", + " \n", + " Parameters:\n", + " u - state at current step \n", + " f - RHS of equation\n", + " dt- time step size\n", + " \n", + " Returns: \n", + " x - array of values at next time step.\n", + " \"\"\"\n", + " \n", + " return u + dt*f(u)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def RK4_DP(u, f, dt):\n", + " \"\"\"Runge Kutta fourth order integration method\n", + " \n", + " Parameters:\n", + " u - state of the system at time t\n", + " f - function for RHS of state equations\n", + " dt - time step\n", + " \n", + " Returns: \n", + " state values at next time step.\n", + " \"\"\"\n", + "\n", + " k1 = dt*f(u)\n", + " k2 = dt*f(u+ k1/2)\n", + " k3 = dt*f(u + k2/2)\n", + " k4 = dt*f(u + k3)\n", + " \n", + " return u + 1/6*(k1 + 2*k2 + 2*k3 + k4)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " this.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width);\n", + " canvas.attr('height', height);\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "#Euler\n", + "x1_dp = numpy.zeros((N,4)) \n", + "x1_dp[0,:] = x_init_dp.copy() #set initial conditions\n", + "for n in range(N-1): #integrate with Euler\n", + " x1_dp[n+1,:] = euler_DP(x1_dp[n,:], f_double_pendulum, dt)\n", + " \n", + "#Runge-Kutta\n", + "x4_dp = numpy.zeros((N,4)) \n", + "x4_dp[0,:] = x_init_dp.copy() #set initial conditions\n", + "for n in range(N-1): #integrate with RK4\n", + " x4_dp[n+1,:] = RK4_DP(x4_dp[n,:], f_double_pendulum, dt)\n", + " \n", + "pyplot.figure(figsize=(10,8));\n", + "pyplot.grid(True);\n", + "pyplot.xlabel(r't', fontsize=18);\n", + "pyplot.ylabel(r'position (radians)', fontsize=18);\n", + "pyplot.title('Double Pendulum Euler');\n", + "pyplot.plot(t, x1_dp[:,0], lw=2, label='Joint 1');\n", + "pyplot.plot(t, x1_dp[:,1], lw=2, label='Joint 2')\n", + "pyplot.legend();\n", + "\n", + "pyplot.figure(figsize=(10,8));\n", + "pyplot.grid(True);\n", + "pyplot.xlabel(r't', fontsize=18);\n", + "pyplot.ylabel(r'position (radians)', fontsize=18);\n", + "pyplot.title('Double Pendulum RK4');\n", + "pyplot.plot(t, x4_dp[:,0], lw=2, label='Joint 1');\n", + "pyplot.plot(t, x4_dp[:,1], lw=2, label='Joint 2')\n", + "pyplot.legend();" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false + }, + "source": [ + "Again, we can see that these integrators cause the system to go unstable. In fact, there is a physical interpretation for what is happening. In the Euler plot, you can see the position of joint 1 oscillating between almost -pi and almost pi (corresponding to straight up in the air), until it finally reaches a value greater than pi. So, once joint 1 'flips around,' the system diverges.\n", + "\n", + "With simple harmonic motion, the symplectic integrator was able to accurately represent the state of the system without adding energy. Let's see if it is successful with the double pendulum. " + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def Verlet_DP(u, dt):\n", + " \"\"\" Verlet integration for integrating from Hamiltonian equations of motion\n", + " for a ddouble pendulum\n", + " Parameters:\n", + " u - state at current time step\n", + " dt - time step size\n", + " Returns: \n", + " state at next time step.\n", + " \"\"\"\n", + " q1 = u[0]\n", + " p1 = u[1]\n", + " q2 = u[2]\n", + " p2 = u[3]\n", + " \n", + " p1_half = p1 - dt/2*(p1*p2/(m1*l1*l2)*sin(q2-q1)+(m1+m2)*g*l1*sin(q1))\n", + " q1_plus = q1 + dt*(-p1_half/I1 - m2*p1_half/(m1*I1) - p2/(m1*l1*l2)*cos(q2-q1))\n", + " p1_plus = p1_half - dt/2*(p1*p2/(m1*l1*l2)*sin(q2-q1_plus)+(m1+m2)*g*l1*sin(q1_plus))\n", + " \n", + " p2_half = p2 - dt/2*( -p1*p2/(m1*l1*l2)*sin(q2-q1) + m2*g*l2*sin(q2))\n", + " q2_plus = q2 + dt*(-p1/(m1*l1*l2)*cos(q2-q1) - p2_half/I2)\n", + " p2_plus = p2_half - dt/2*( -p1*p2/(m1*l1*l2)*sin(q2_plus-q1) + m2*g*l2*sin(q2_plus))\n", + " \n", + " return numpy.array([q1_plus, p1_plus, q2_plus, p2_plus]) " + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " this.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width);\n", + " canvas.attr('height', height);\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "Verlet_energy = numpy.zeros_like(t)\n", + "Verlet_vel1 = xv_dp[:,1]/m1\n", + "Verlet_vel2 = xv_dp[:,3]/m2\n", + "RK4_energy = numpy.zeros_like(t)\n", + "\n", + "\n", + "u = numpy.array([xv_dp[:,0], Verlet_vel1, xv_dp[:,2], Verlet_vel2])\n", + "#print(u[0,:])\n", + "for i in range(N):\n", + " Verlet_energy[i] = get_Energy(u[:,i])\n", + " RK4_energy[i] = get_Energy(x4_dp[i,:])\n", + "\n", + " \n", + "\n", + "pyplot.figure(figsize=(10,8));\n", + "pyplot.grid(True);\n", + "\n", + "pyplot.title('Energy in System');\n", + "ax1 = pyplot.gca()\n", + "pyplot.plot(t, Verlet_energy, lw=2, label='Verlet', color='b');\n", + "pyplot.xlabel('time (s)')\n", + "pyplot.ylabel('Energy')\n", + "pyplot.plot(t, RK4_energy, 'r', lw=2, label='RK4')\n", + "pyplot.legend();\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "So, we are very far from perfect. But, at least with Verlet the energy in the system is stable, even if it isn't constant! " + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "### Investigating system behavior\n", + "\n", + "In order to investigate system behavior for damped harmonic motion, we used a phase portrait where we plotted position vs. velocity. In this system, we have 2 positions and 2 velocities, so a phase plot won't accurately capture the behavior of the system. Instead we are going to use something called a Poincare section (https://en.wikipedia.org/wiki/Poincar%C3%A9_map). The way these work is to capture the system behavior at a specific, repetetive point. For example, each time that $\\theta_2 = 0$ or any multiple of 2*pi (mass 2 is hanging straight down) we can capture the angular position and velocity of mass 2. This way, we are taking a 'snapshot' of the system, and plotting the state each time. We can now investigate system behavior on a simple 2D plot, similar to the phase portrait." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def poincare(state, sf, val, s1, s2):\n", + " '''Creates a poincare section, capturing system state each time q[state]=val\n", + " \n", + " Parameters:\n", + " state - state vector\n", + " sf - the row in q to fix\n", + " val - the value to capture state at \n", + " s1 - state to capture 1\n", + " s2 - state to capture 2\n", + " \n", + " Returns:\n", + " v - array with information for poincare section\n", + " '''\n", + " v = numpy.zeros((2,1))\n", + " for i in range(len(state[:,sf])-1):\n", + " if (state[i,sf]>val):\n", + " if (state[i+1,sf] < val):\n", + " v = numpy.append(v, [[state[i, s1]],[state[i, s2]]], axis=1)\n", + " \n", + " if (state[i,sf] < val):\n", + " if (state[i+1,sf] > val):\n", + " v = numpy.append(v, [[state[i, s1]],[state[i, s2]]], axis=1)\n", + " \n", + " return v" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " this.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width);\n", + " canvas.attr('height', height);\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " this.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width);\n", + " canvas.attr('height', height);\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "v = poincare(q4d_dp, 0, 0, 1, 3) #use states from RK4 integration\n", + "pyplot.figure(figsize=(10,8));\n", + "pyplot.grid(True);\n", + "pyplot.xlabel(r'$\\theta_2$', fontsize=18);\n", + "pyplot.ylabel(r'$\\dot{\\theta_2}$', fontsize=18);\n", + "pyplot.title('Poincare Section');\n", + "pyplot.plot(v[0,:], v[1,:], 'bo', lw=2);\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Other than a few outliers, the plot converges seems to converge. So, this system is stable! This would not be the case for our undamped double pendulum, as integrated by RK4. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Sources: \n", + "\n", + "[1] http://scienceworld.wolfram.com/physics/DoublePendulum.html\n", + "\n", + "[2] http://www.phy.uct.ac.za/courses/opencontent/phylab2/worksheet9_09.pdf\n", + "\n", + "[3] http://www.phys.lsu.edu/faculty/gonzalez/Teaching/Phys7221/DoublePendulum.pdf\n", + "\n", + "[4] http://www.iontrap.wabash.edu/adlab/papers/F2011_foster_groninger_tang_chaos.pdf\n", + "\n", + "[5] https://math.berkeley.edu/~alanw/242papers99/markiewicz.pdf\n", + "\n", + "[6] http://www.unige.ch/~hairer/poly_geoint/week2.pdf\n", + "\n", + "[7] http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?1994CeMDA..60..409T&defaultprint=YES&filetype=.pdf\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# This cell loads the style of the notebook, which is modified from the \n", + "# Numerical Methods in Python Course: http://openedx.seas.gwu.edu/courses/GW/MAE6286/2014_fall/about\n", + "\n", + "from IPython.core.display import HTML\n", + "css_file = './numericalmoocstyle.css'\n", + "HTML(open(css_file, \"r\").read())" + ] + } + ], + "metadata": { + "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.5.0" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/randy_schur/Integration_Techniques_for_Mechanical_Systems.ipynb b/randy_schur/Integration_Techniques_for_Mechanical_Systems.ipynb new file mode 100644 index 0000000..0ee3e41 --- /dev/null +++ b/randy_schur/Integration_Techniques_for_Mechanical_Systems.ipynb @@ -0,0 +1,7546 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Integration of Mechanical Systems\n", + "\n", + "In this module, we will be exploring numerical integrators and how they perform on different mechanical systems." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Simple Harmonic Oscillator\n", + "\n", + "The simple harmonic oscillator is a simple physical system described by the second-order differential equation below. Despite its simplicity, this is a system that shows up in a similar form in many different fields of engineering. The same equation describes the behavior of a spring-mass-damper system, RLC circuit, or the motion of a (small-angle) pendulum.\n", + "\n", + "![Image](figures/spring-mass-damper.png)\n", + "\n", + "##### Figure 1. Diagram of Spring-Mass-Damper system (from https://en.wikipedia.org/wiki/Harmonic_oscillator)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The differential equation describing the motion of a simple harmonic oscillator is:\n", + "\n", + "$$\\begin{equation}\n", + "0 = m\\ddot{x} + kx\n", + "\\end{equation}$$\n", + "\n", + "Since this is a second order system, we will need two states to accurately describe the system. We can use position $x$ and velocity $\\dot{x}$. This is called the state vector $q$, as:\n", + "\n", + "$$\\begin{eqnarray}\n", + "\\vec{q} = \\begin{bmatrix} x\\\\ \\dot{x} \\end{bmatrix}\n", + "\\end{eqnarray}$$\n", + "\n", + "This system can then be broken into two first order differential equations:\n", + "\n", + "$$\\begin{eqnarray}\n", + "\\vec{\\dot{q}} = f(\\vec{q}) = \\begin{bmatrix} q_2\\\\ -\\frac{k}{m}q_1 \\end{bmatrix}\n", + "\\end{eqnarray}$$\n", + "\n", + "To look at the behavior of this system, we can start with some initial conditions for $q$, then use a numerical integration method such as Euler's method to integrate forward in time. This is equivalent to pulling (or pushing) on the spring, then letting go at t=0." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Equations of Motion from Lagrangian\n", + "\n", + "You are encouraged to follow this derivation ON PAPER on your own. The equations of motion are given above for a simple harmonic oscillator, but what if we didn't have this equation ahead of time? The following is a method that will work in nearly any situation to find equations of motion, often working well where a force analysis becomes difficult. This is especially good practice for anyone studying mechanics. To derive the equations of motion, we will start with the Lagrangian. This is defined by:\n", + "\n", + "$$\\begin{equation}\n", + "L = T- U\n", + "\\end{equation}$$\n", + "\n", + "where $T$ is the total Kinetic Energy in the system, and $U$ is the total Potential Energy. Potential energy in this system comes from the stretching or compressing of the spring, and is calculated as:\n", + "\n", + "$$\\begin{eqnarray*}\n", + "U &=& \\frac{1}{2}kx^2\\\\\n", + "\\end{eqnarray*}$$\n", + "\n", + "Kinetic Energy in this system is from the motion of the mass, and is calculated as:\n", + "\n", + "$$\\begin{eqnarray*}\n", + "T &=& \\frac{1}{2}mv^2 \\\\\n", + "&=& \\frac{1}{2}m\\dot{x}^2 \n", + "\\end{eqnarray*}$$\n", + "\n", + "\n", + "So the Lagrangian quantity becomes:\n", + "$$\\begin{equation}\n", + "L = \\frac{1}{2}m\\dot{x}^2 - \\frac{1}{2}kx^2\\\\\n", + "\\end{equation}$$\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The equations of motion are then found using the Lagrange Equation:\n", + "\n", + "$$\\begin{equation}\n", + "\\frac{d}{dt} \\left( \\frac{\\partial L}{\\partial \\dot{x_i}}\\right) - \\frac{\\partial L}{\\partial x_i} = 0\n", + "\\end{equation}$$\n", + "\n", + "Since we have only one generalized coordinate (position), we need only one equation of motion. We'll calculate this in steps:\n", + "\n", + "$$\\begin{eqnarray*}\n", + "\\frac{\\partial L}{\\partial \\dot{x}} &=& m\\dot{x} \\\\\n", + "\\frac{d}{dt} \\left( \\frac{\\partial L}{\\partial \\dot{x}}\\right) &=& m\\ddot{x}\\\\\n", + "-\\frac{\\partial L}{\\partial x} &=& + kx\n", + "\\end{eqnarray*}$$\n", + "\n", + "We can add together the second and third equations to find the first equation of motion. Some of these terms will cancel out, and we are left with a final equation of motion:\n", + "\n", + "$$\\begin{equation}\n", + "0 = m\\ddot{x} + kx\n", + "\\end{equation}$$\n", + "\n", + "This is the same equation we have before! So, we can be relatively confident that this derivation is correct." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Integration Methods\n", + "\n", + "### Euler's method\n", + "\n", + "Euler's method is discussed in lesson 1 of the Numerical Methods course. The method goes as follows:\n", + "$$x_i^{n+1} = x_i^n + h\\dot{x}_i^n$$\n", + "\n", + "where $h$ is the timestep. This isn't a complicated scheme, so for details see the code below." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Runge-Kutta Integration (RK4)\n", + "\n", + "Depending on the application, Euler's method may provide enough accuracy with a small timestep, especially a simple harmonic oscillator. However, Euler is a first-order method and we can do better. One option is a Runge-Kutta scheme. This is a popular numerical integration method which is easily extended to higher orders. Here we will use a fourth order Runge-Kutta, also called RK4:\n", + "\n", + "$$\\begin{equation}\n", + "q_i^{n+1} = q_i^n + \\frac{h}{6}\\left(k_1+2k_2+2k_3+k_4\\right)\n", + "\\end{equation}$$\n", + "where\n", + "$$\\begin{eqnarray}\n", + "k_1 &=& f(t_n, x_n)\\\\\n", + "k_2 &=& f(t_n+\\frac{h}{2}, x_n+\\frac{h}{2}k_1)\\\\\n", + "k_3 &=& f(t_n+\\frac{h}{2}, x_n+\\frac{h}{2}k_2)\\\\\n", + "k_4 &=& f(t_n+h, x_n+hk_3)\\\\\n", + "t_{n+1} &=& t_n + h\n", + "\\end{eqnarray}$$\n", + "\n", + "For more information on the Runge-Kutta method, including the derivation and explanation of coefficients, see [9].\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Symplectic Integrators\n", + "\n", + "Symplectic integrators are similar to the methods discussed above, but they use equations of motion derived from Hamiltonian mechanics. According to [5], syplectic integrators preserve the conserved Hamiltonian quantities. In practical terms, this works out to mean the methods reflect conservation of energy, down to a truncation error.\n", + "\n", + "The Hamiltonian equations of motion can be derived from the total energy in the system, much like the Lagrangian. The coordinates we use are position ($q$) and momentum ($p = m\\dot x$). The Hamiltonian is:\n", + "$$ H = T+V$$\n", + "where $T$ is the kinetic energy in the system, and $V$ is the potential energy in the system. So, we get:\n", + "$$\\begin{eqnarray*}\n", + "T &=& \\frac{1}{2}m\\dot{x}^2 = \\frac{p^2}{2m}\\\\\n", + "V &=& \\frac{1}{2}k{x}^2 = \\frac{kq^2}{2} \n", + "\\end{eqnarray*}$$\n", + "\n", + "and \n", + "\n", + "$$H = \\frac{p^2}{2m} + \\frac{kq^2}{2}$$\n", + "\n", + "\n", + "\n", + "The Hamilton's equations are:\n", + "$$\\begin{eqnarray*}\n", + "\\frac{dq}{dt} &=& -\\frac{\\partial H}{\\partial p}\\\\\n", + "\\frac{dp}{dt} &=& \\frac{\\partial H}{\\partial q}\n", + "\\end{eqnarray*}$$\n", + "\n", + "So our equations of motion become:\n", + "$$\\begin{equation*}\n", + "\\begin{bmatrix} \\dot q \\\\ \\dot p \\end{bmatrix} = \\begin{bmatrix} -\\frac{\\partial H}{\\partial p}\\\\ \\frac{\\partial H}{\\partial q} \\end{bmatrix} = \\begin{bmatrix} -\\frac{p}{m} \\\\ kq\\end{bmatrix}\n", + "\\end{equation*}$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's now take a look at symplectic integration methods [6]. We'll start with a first order method - symplectic Euler's method. This uses the equations:\n", + "\n", + "$$\\begin{eqnarray*}\n", + "p_{n+1} = p_n - h \\frac{\\partial H}{\\partial q} \\Bigr|_{p_{n+1},q_n}\\\\\n", + "q_{n+1} = q_n + h\\frac{\\partial H}{\\partial p} \\Bigr|_{p_{n+1},q_n}\n", + "\\end{eqnarray*}$$\n", + "\n", + "We can then plug in the partial derivatives and integrate the system. This system will do a better job of conserving energy than Euler's method (we'll investigate this later), $\\textit{but only to a truncation error.}$ Again, we can do better than first-order. Let's try a second-order symplectic scheme, also called Verlet integration. We'll have the same equations of motion, but this time a different set of integration equations [6]:\n", + "\n", + "$$\\begin{eqnarray*}\n", + "p_{n+1/2} &=& p_n - \\frac{h}{2} \\frac{\\partial H}{\\partial q} \\Bigr|_{p_{n+1/2},q_n}\\\\\n", + "q_{n+1} &=& q_n + \\frac{h}{2} \\left(\\frac{\\partial H}{\\partial p} \\Bigr|_{p_{n+1/2},q_n} + \\frac{\\partial H}{\\partial p} \\Bigr|_{p_{n+1/2},q_{n+1/2}} \\right) \\\\\n", + "p_{n+1} &=& p_{n+1/2} - \\frac{h}{2} \\frac{\\partial H}{\\partial q}\\Bigr|_{p_{n+1/2},q_{n+1/2}}\n", + "\\end{eqnarray*}$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Analytical Solution\n", + "The simple harmonic oscillator is a nice example for numerical integration, because it has an analytical form which we can compare to. The analytical solution is\n", + "\n", + "$$\\begin{equation}\n", + "x(t) = Acos(\\omega_n t + \\phi)\n", + "\\end{equation}$$\n", + "\n", + "where $\\omega_n$ is the natural frequency of the system given by $\\omega_n = \\sqrt{\\frac{k}{m}}$ . The amplitude, $A$, and the phase, $\\phi$, of oscillation are determined from the initial conditions. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Let's Get to Programming" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import numpy\n", + "from scipy.linalg import solve\n", + "from numpy.linalg import det\n", + "from math import pi, cos, sin, sqrt\n", + "\n", + "from matplotlib import pyplot\n", + "from matplotlib.pyplot import quiver\n", + "%matplotlib notebook\n", + "from matplotlib import rcParams, cm\n", + "rcParams['font.family'] = 'serif'\n", + "rcParams['font.size'] = 16" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "#Set up system parameters:\n", + "m = 12\n", + "k = 150\n", + "\n", + "T = 25; #[seconds]\n", + "dt = .02; #time step size\n", + "N = int(T/dt)+1 #number of time steps\n", + "t = numpy.linspace(0.0, T, N) #array of time values" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def euler(u, f, dt):\n", + " \"\"\" Euler's method for integrating a system of differential equations.\n", + " \n", + " Parameters:\n", + " u - state at current step \n", + " f - RHS of equation\n", + " dt- time step size\n", + " \n", + " Returns: \n", + " x - array of values at next time step.\n", + " \"\"\"\n", + " \n", + " return u + dt*f(u) " + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def RK2(u, f, dt):\n", + " \"\"\"Runge Kutta second order integration method\n", + " \n", + " Parameters:\n", + " u - state of the system at time t\n", + " f - function for RHS of state equations\n", + " dt - time step\n", + " \n", + " Returns: \n", + " state values at next time step.\n", + " \"\"\"\n", + " k1 = u + 0.5*dt*f(u)\n", + " k2 = u + dt*f(k1)\n", + " return k2" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def RK4(u, f, dt):\n", + " \"\"\"Runge Kutta fourth order integration method\n", + " \n", + " Parameters:\n", + " u - state of the system at time t\n", + " f - function for RHS of state equations\n", + " dt - time step\n", + " \n", + " Returns: \n", + " state values at next time step.\n", + " \"\"\"\n", + " k1 = dt*f(u)\n", + " k2 = dt*f(u+ k1/2)\n", + " k3 = dt*f(u + k2/2)\n", + " k4 = dt*f(u + k3)\n", + " \n", + " return u + 1/6*(k1 + 2*k2 + 2*k3 + k4)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def f_harmonic_oscillator(u):\n", + " \"\"\"Returns RHS of harmonic oscillator EOM\n", + " \n", + " Parameters:\n", + " u - initial state\n", + " \n", + " Returns:\n", + " RHS of harmonic oscillator eqn.\n", + " \n", + " \"\"\"\n", + " pos = u[0]\n", + " vel = u[1]\n", + " \n", + " return numpy.array([vel, -k/m*pos])" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def verlet_SHM(u, dt):\n", + " \"\"\" Verlet integration for integrating from Hamiltonian equations of motion\n", + " for a simple harmonic oscillator\n", + " Parameters:\n", + " u - state at current time step\n", + " dt - time step size\n", + " Returns: \n", + " state at next time step.\n", + " \"\"\"\n", + " pos = u[0]\n", + " mom = u[1]\n", + " \n", + " q_half = pos - dt/2*mom/m\n", + " #p_half = mom + dt/2*(k*pos)\n", + " \n", + " p = mom + dt/2*(2*k*q_half)\n", + " q = q_half - dt/2*p/m\n", + " \n", + " return numpy.array([q, p]) " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "#Initial Conditions\n", + "x0 = 0.5 #[m]\n", + "xdot0 = 0 #[m/s]\n", + "\n", + "#Equivalent initial conditions for the Hamiltonian system\n", + "q0 = x0\n", + "p0 = m*xdot0\n", + "\n", + "x_init = numpy.array([x0, xdot0])\n", + "x_init_H = numpy.array([q0, p0])" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " this.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width);\n", + " canvas.attr('height', height);\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pyplot.figure(figsize=(10,6))\n", + "pyplot.tick_params(axis='both', labelsize=14)\n", + "pyplot.grid(True)\n", + "pyplot.xlabel('$\\Delta t$', fontsize=16)\n", + "pyplot.ylabel('Error', fontsize=16)\n", + "pyplot.loglog(dt_values, Euler_error_values, 'bo-', label='Euler')\n", + "pyplot.loglog(dt_values, RK2_error_values, 'ro-', label='RK2')\n", + "pyplot.loglog(dt_values, RK4_error_values, 'co-', label='RK4')\n", + "pyplot.loglog(dt_values, Verlet_error_values, 'go-', label='Symplectic')\n", + "pyplot.axis('equal')\n", + "pyplot.title('Error in Simple Harmonic Motion Integration')\n", + "pyplot.legend();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Is this what you would expect to see based on the plots of position above? What are the sources of this error? We have the usual suspects: truncation error, discretization error, possible precision errors. Don't believe that these are significant? To examine precision errors, take a look at the following screenshot, which shows the plots above, built with identical code on two different computers. The responses on the right are from a 64-bit computer, and the responses on the left are from a 32-bit computer. Look closely at the phase shift in the system response. That's a huge difference! Even the precision of the machine can have a significant impact on numerical integration.\n", + "\n", + "![Image](figures/combined_response.png)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Despite the numerical errors, why is one second order scheme (Verlet) outperforming another (RK2) by so much? What is going on? We can check the order of convergence to make sure that these schemes are actually 2nd order." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/mint/miniconda3/lib/python3.5/site-packages/ipykernel/__main__.py:24: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n" + ] + } + ], + "source": [ + "def get_diffgrid(u_current, u_fine, dt):\n", + " \"\"\"Returns the difference between one grid and the fine one using L-1 norm.\n", + " \n", + " Parameters\n", + " ----------\n", + " u_current : array of float\n", + " solution on the current grid.\n", + " u_finest : array of float\n", + " solution on the fine grid.\n", + " dt : float\n", + " time-increment on the current grid.\n", + " \n", + " Returns\n", + " -------\n", + " diffgrid : float\n", + " difference computed in the L-1 norm.\n", + " \"\"\"\n", + " \n", + " N_current = len(u_current[:,0])\n", + " N_fine = len(u_fine[:,0])\n", + " \n", + " grid_size_ratio = numpy.ceil(N_fine/N_current)\n", + " \n", + " diffgrid = dt * numpy.sum( numpy.abs(\\\n", + " u_current[:,0]- u_fine[::grid_size_ratio,0])) \n", + " \n", + " return diffgrid\n", + "\n", + "# use a for-loop to compute the solution on different grids - this code is modified from Numerical MOOC lesson 1 module 4\n", + "dt_values = numpy.array([0.1, 0.05, 0.01, 0.005, 0.001])\n", + "u1_values = numpy.empty_like(dt_values, dtype=numpy.ndarray)\n", + "u2_values = u1_values.copy()\n", + "u4_values = u1_values.copy()\n", + "uv_values = u1_values.copy()\n", + "\n", + "#At each value of dt, integrate the system using 4 different methods from t=0 to T\n", + "for i, dt in enumerate(dt_values):\n", + " N = int(T/dt)+1 \n", + " t = numpy.linspace(0.0, T, N)\n", + " u1 = numpy.empty((N,2))\n", + " u1[0] = numpy.array([x0, xdot0])\n", + " u2 = u1.copy()\n", + " u4 = u1.copy()\n", + " uv = u1.copy()\n", + " for n in range(N-1):\n", + " u1[n+1] = euler(u1[n], f_harmonic_oscillator, dt)\n", + " u2[n+1] = RK2(u2[n], f_harmonic_oscillator, dt)\n", + " u4[n+1] = RK4(u4[n], f_harmonic_oscillator, dt)\n", + " uv[n+1] = verlet_SHM(uv[n], dt)\n", + " # store the value of u related to one grid\n", + " u1_values[i] = u1\n", + " u2_values[i] = u2\n", + " u4_values[i] = u4\n", + " uv_values[i] = uv\n", + " \n", + "# compute diffgrid\n", + "diffgrid1 = numpy.empty_like(dt_values)\n", + "diffgrid2 = diffgrid1.copy()\n", + "diffgrid4 = diffgrid1.copy()\n", + "diffgridv = diffgrid1.copy()\n", + "for i, dt in enumerate(dt_values):\n", + " diffgrid1[i] = get_diffgrid(u1_values[i], u1_values[-1], dt)\n", + " diffgrid2[i] = get_diffgrid(u2_values[i], u2_values[-1], dt)\n", + " diffgrid4[i] = get_diffgrid(u4_values[i], u4_values[-1], dt)\n", + " diffgridv[i] = get_diffgrid(uv_values[i], uv_values[-1], dt)" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " this.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width);\n", + " canvas.attr('height', height);\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " this.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width);\n", + " canvas.attr('height', height);\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "#Euler\n", + "x1_d = numpy.zeros((N,2)) \n", + "x1_d[0] = x_init.copy() #set initial conditions\n", + "#Runge-Kutta Second Order\n", + "x2_d = x1_d.copy()\n", + "#Runge-Kutta Fourth Order\n", + "x4_d= x1_d.copy()\n", + "#Symplectic Second Order\n", + "xv_d = x1_d.copy()\n", + "\n", + "for n in range(N-1): #integrate with Euler\n", + " x1_d[n+1] = euler(x1_d[n], f_damped_HM, dt)\n", + " x2_d[n+1] = RK2(x2_d[n], f_damped_HM, dt)\n", + " x4_d[n+1] = RK4(x4_d[n], f_damped_HM, dt)\n", + " xv_d[n+1] = Verlet_damped_HM(xv_d[n], dt)\n", + " \n", + "#Analytical Solution \n", + "A = x0 #with no forcing, the max amplitude is equal to initial amplitude\n", + "phi = 0\n", + "x_d_analytical = A*numpy.exp(-zeta*wn*t)*numpy.cos(wd*t + phi)\n", + "vel_d_analytical = -A*numpy.exp(-zeta*wn*t)*wd*numpy.sin(wd*t+phi)\n", + " \n", + "pyplot.figure(figsize=(10,8));\n", + "pyplot.grid(True);\n", + "pyplot.xlabel(r't', fontsize=18);\n", + "pyplot.ylabel(r'position (meters)', fontsize=18);\n", + "pyplot.title('Harmonic oscillator position');\n", + "pyplot.plot(t, x1_d[:,0], lw=2, label='Euler');\n", + "pyplot.plot(t, x2_d[:,0], 'r-', lw=2, label='RK2');\n", + "pyplot.plot(t, x4_d[:,0], 'c-', lw=2, label='RK4');\n", + "pyplot.plot(t, xv_d[:,0], 'g', lw=2, label='Verlet');\n", + "pyplot.plot(t, x_d_analytical, 'k--', lw=2, label='analytical');\n", + "pyplot.legend();\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*Note: In the figure above, the energy in the system according to Verlet, RK2, and RK4 is so close that the lines are appearing on top of each other. If you zoom in really closely you can see the differences. \n", + "\n", + "How does each integrator perform compared to the case of simple harmonic motion? Just by looking at the system response, it seems like the error is actually better here. Why is that? It is important to note here that the integrators perform differently for different system models. This is part of why we study numerical methods! There isn't a single answer to which scheme to use, it is dependent on the characteristics of both the scheme and the system.\n", + "\n", + "The effect of the damper is to remove energy from the system. So, in this case the damping term is helping to remove some energy from the system that the integration adds. As time goes to infinity, these solutions will converge, and the Euler and Runge-Kutta methods are stable. We can again look at the energy in the system according to each method. \n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " this.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width);\n", + " canvas.attr('height', height);\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "#create a phase plot\n", + "NP = 15\n", + "pos = numpy.linspace(-x0, x0, NP)\n", + "vel = numpy.linspace(-2*x0, 2*x0, NP)\n", + "\n", + "X1, X2 = numpy.meshgrid(pos, vel)\n", + "u = numpy.zeros_like(X1)\n", + "v = numpy.zeros_like(X2)\n", + "\n", + "for i in range(NP):\n", + " for j in range(NP):\n", + " x = X1[i, j]\n", + " y = X2[i, j]\n", + " yprime = f_damped_HM([x, y])\n", + " u[i,j] = yprime[0]\n", + " v[i,j] = yprime[1]\n", + "\n", + "\n", + "pyplot.figure(figsize=(10,8));\n", + "pyplot.quiver(X1, X2, u, v, pivot='mid', color='k', label='direction field');\n", + "pyplot.plot(x_d_analytical, vel_d_analytical, 'k-', label='analytical');\n", + "pyplot.plot(x1_d[:,0], x1_d[:,1],'b-', label='Euler');\n", + "pyplot.plot(x2_d[:,0], x2_d[:,1], 'r-', label='RK2');\n", + "pyplot.plot(x4_d[:,0], x4_d[:,1], 'c-', label='RK4');\n", + "pyplot.plot(xv_d[:,0], -Verlet_vel_d, 'g-', label='Verlet');\n", + "pyplot.xlabel('Postion ($x_1$)')\n", + "pyplot.ylabel('Velocity ($x_2$)')\n", + "pyplot.legend();\n", + "pyplot.title('Phase portrait for damped harmonic oscillator');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this case, it appears that Verlet, RK2, and RK4 give almost identical results, which are close but not exactly the same as the analytical solution. So for the damped harmonic oscillator, there is less of an advantage to using a more complex integration scheme. Try modifying the code above to make a phase portrait for the simple harmonic oscillator. We should see a circle (or an ellipse) form because the system behavior repeats itself. This is called a limit cycle. Which integrators create this behavior?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Conclusion\n", + "\n", + "In this notebook we discussed four different numerical integrators and applied these to a mechanical system. Even for the relatively straightforward case of harmonic motion, we found that different integrators can have very different behavior. In order to choose the best method (not too complex, but still accurate) we investigated how each method affects conservation of energy. We also found that not every situation is the same - the integrators performed differently for damped and undamped systems.\n", + "\n", + "In the next notebook, we will take these numerical integrators and apply them to the more complex double pendulum system. We'll use some of the same tools to analyze the behavior of each numerical integrator. However, there is no closed form solution for the double pendulum, so we'll use what we learned in this notebook to try and determine the effectiveness of each method on the more complex system." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false + }, + "source": [ + "Sources and additional resources: \n", + "\n", + "[1] http://scienceworld.wolfram.com/physics/DoublePendulum.html\n", + "\n", + "[2] http://www.phy.uct.ac.za/courses/opencontent/phylab2/worksheet9_09.pdf\n", + "\n", + "[3] http://www.phys.lsu.edu/faculty/gonzalez/Teaching/Phys7221/DoublePendulum.pdf\n", + "\n", + "[4] http://www.iontrap.wabash.edu/adlab/papers/F2011_foster_groninger_tang_chaos.pdf\n", + "\n", + "[5] https://math.berkeley.edu/~alanw/242papers99/markiewicz.pdf\n", + "\n", + "[6] http://www.unige.ch/~hairer/poly_geoint/week2.pdf\n", + "\n", + "[7] http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?1994CeMDA..60..409T&defaultprint=YES&filetype=.pdf\n", + "\n", + "[8] http://www.ujfi.fei.stuba.sk/fyzika-ballo/Runge-Kutta.pdf\n", + "\n", + "[9] http://web.mit.edu/10.001/Web/Course_Notes/Differential_Equations_Notes/node5.html" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# This cell loads the style of the notebook, which is modified from the \n", + "# Numerical Methods in Python Course: http://openedx.seas.gwu.edu/courses/GW/MAE6286/2014_fall/about\n", + "\n", + "from IPython.core.display import HTML\n", + "css_file = './numericalmoocstyle.css'\n", + "HTML(open(css_file, \"r\").read())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "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.5.0" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/randy_schur/figures/Error_64bit.png b/randy_schur/figures/Error_64bit.png new file mode 100644 index 0000000..2cb22b4 Binary files /dev/null and b/randy_schur/figures/Error_64bit.png differ diff --git a/randy_schur/figures/POS_computer.png b/randy_schur/figures/POS_computer.png new file mode 100644 index 0000000..65c27cc Binary files /dev/null and b/randy_schur/figures/POS_computer.png differ diff --git a/randy_schur/figures/POS_computer2.png b/randy_schur/figures/POS_computer2.png new file mode 100644 index 0000000..0f25ac8 Binary files /dev/null and b/randy_schur/figures/POS_computer2.png differ diff --git a/randy_schur/figures/combined_response.png b/randy_schur/figures/combined_response.png new file mode 100644 index 0000000..ceb716c Binary files /dev/null and b/randy_schur/figures/combined_response.png differ diff --git a/randy_schur/figures/diagram.png b/randy_schur/figures/diagram.png new file mode 100644 index 0000000..a613042 Binary files /dev/null and b/randy_schur/figures/diagram.png differ diff --git a/randy_schur/figures/error_combined.png b/randy_schur/figures/error_combined.png new file mode 100644 index 0000000..c78925f Binary files /dev/null and b/randy_schur/figures/error_combined.png differ diff --git a/randy_schur/figures/response_64bit.png b/randy_schur/figures/response_64bit.png new file mode 100644 index 0000000..3a6e320 Binary files /dev/null and b/randy_schur/figures/response_64bit.png differ diff --git a/randy_schur/figures/spring-mass-damper.png b/randy_schur/figures/spring-mass-damper.png new file mode 100644 index 0000000..6d473b7 Binary files /dev/null and b/randy_schur/figures/spring-mass-damper.png differ diff --git a/randy_schur/numericalmoocstyle.css b/randy_schur/numericalmoocstyle.css new file mode 100644 index 0000000..e6a8035 --- /dev/null +++ b/randy_schur/numericalmoocstyle.css @@ -0,0 +1,142 @@ + + + + + + +