{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# odeint examples\n", "\n", "Code for *Modeling and Simulation in Python*\n", "\n", "by Allen B. Downey, available from http://greenteapress.com\n", "\n", "Copyright 2017 Allen B. Downey\n", "\n", "MIT License: https://opensource.org/licenses/MIT" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import print_function, division\n", "\n", "#%matplotlib notebook\n", "%matplotlib inline\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "\n", "import numpy as np\n", "from scipy.integrate import odeint\n", "\n", "from pint import UnitRegistry\n", "UNITS = UnitRegistry()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Basic example using odeint without units\n", "\n", "Here is a slope function for a spring mass system where the spring force depends on positive and the friction force depends on velocity." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def slope_func(Y, t):\n", " y1, y2 = Y\n", " y1p = y2\n", " y2p = -(y1 + y2)\n", " return y1p, y2p" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Starting with initial position 1 and velocity 0." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[1, 0]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y_init = [1, 0]\n", "y_init" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we evaluate the slope function with the starting conditions, we get velocity toward 0." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "(0, -1)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "slope_func(y_init, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's the sequence of times we would like to integrate over." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "numpy.ndarray" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ts = np.arange(0, 15.0, 0.1)\n", "type(ts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here's how we call `odeint`:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "asol = odeint(slope_func, y_init, ts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result is a matrix, but we can assign the columns to variables." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "position, velocity = asol.transpose()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And plot position as a function of time." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG9ZJREFUeJzt3Xm0VOWZ7/Hvc0BAQZEhIoMgSBSDoiGKKIglJgu0VbLs\n1oa4YjRta27EOKRvNLnJAnt1J+pq1zVGTK5T2nSixA4m2okDJnqMRlRUZDwIijkcBgEVjaggw3P/\neOtAcThz7VPvrtq/z1p7Ve06m70fpme/9bzDNndHRESypSp2ACIiUnpK/iIiGaTkLyKSQUr+IiIZ\npOQvIpJBSv4iIhmUSPI3s3vMbIOZLWrmmNvMbKWZvWZmxydxXRERaZ+kWv4/ByY19UMzOxM4wt0/\nC1wO/Cyh64qISDskkvzd/TlgczOHTAF+kT/2RaCnmfVL4toiItJ2par5DwTqCvbX5j8TEZEI1OEr\nIpJBnUt0nbXAYQX7g/Kf7cPMtNiQiEgbubu15fgkW/6W3xrzCHARgJmNBd539w1NnWj+fKdvX2fR\nIsc9fduMGTOix6A4FafiVJz1W3sk0vI3s/uBHNDHzFYDM4AugLv7ne7+qJmdZWZvAB8BlzR3vhNO\ngFtvhSlT4OWXoXfvJKIUEZF6iSR/d/9KK46Z3pZzXnghPPMM3Hgj3Hxz+2MTEZF9pbrDd8YMuPtu\nWL8+diR7y+VysUNoFcWZLMWZLMUZl7W3XtRRzMwLY7r2WtixA267LWJQIiIpZmZ4Gzt8U5/8N26E\nESPgtddg8OCIgYmIpFR7kn+qyz4AhxwCF18Md9wROxIRkcqR+pY/QE0NTJwIq1fDfvtFCkxEJKUq\nsuUPcPTRMGwYPPpo7EhERCpDWSR/gEsvDSN/RESkeGVR9gHYsgUOOwyWLIGBWhJORGS3ii37APTo\nAeefD7/8ZexIRETKX9kkf4ALLoA5c2JHISJS/sqm7AOwfTsceigsWKAx/yIi9Sq67ANhmOe558Lv\nfhc7EhGR8lZWyR/gvPPgoYdiRyEiUt7KquwDsHVrKP2sWBFm/4qIZF3Fl30AunWDyZPh4YdjRyIi\nUr7KLvlDqPv/4Q+xoxARKV9lV/YB2LQJhg8Pr126lCgwEZGUykTZB+Azn4GjjoLnn48diYhIeSrL\n5A+h7v/447GjEBEpT2Wb/CdNgieeiB2FiEh5KsuaP4RHO37mM7BsGfTvX4LARERSKjM1f4DOneGL\nX4S5c2NHIiJSfso2+YNKPyIi7VW2ZR+AVatg/HhYuxasTV94REQqR6bKPgBDh4bF3lasiB2JiEh5\nKevkbwa5HDz9dOxIRETKS1knf4DTT4fq6thRiIiUl7Ku+QPU1sJJJ8H69ar7i0g2Rav5m9lkM1tu\nZivM7LpGfn6QmT1iZq+Z2WIzuziJ6wIMGQIHHAA1NUmdUUSk8hWd/M2sCrgdmASMBKaZ2YgGh10B\nLHX344HTgVvMrHOx166nur+ISNsk0fIfA6x091p33w7MBqY0OMaBA/PvDwTedfcdCVwbCHV/JX8R\nkdZLIvkPBOoK9tfkPyt0O/A5M1sHLASuSuC6u02YAM89BynrvhARSa1SjfaZBCxw9wHA54FZZtYj\nqZMPHhzG+7/5ZlJnFBGpbEnU3dcCgwv2B+U/K3QJ8CMAd3/TzN4CRgAvN3bCmTNn7n6fy+XI5XLN\nBmAWZvo+91x4yIuISCWrrq6musgx7kUP9TSzTsDrwBnAeuAlYJq71xQcMwvY6O43mFk/QtI/zt3f\na+R8bRrqWW/WLFiwAO6+u52/ERGRMhVlqKe77wSmA3OBpcBsd68xs8vN7LL8Yf8GnGJmi4Ange80\nlviLUd/yFxGRlpX9JK96O3dCnz6wcmVY519EJCsyt7BboU6d4OST4S9/iR2JiEj6VUzyB5V+RERa\nS8lfRCSDKqbmD/Dxx6He/9570LVrwoGJiKRUpmv+EBZ4O+qoMORTRESaVlHJH2DsWHjhhdhRiIik\nW0Um/xdfjB2FiEi6VWTyV8tfRKR5FZf8P/tZ+OADePvt2JGIiKRXxSV/s/BYR5V+RESaVnHJH1T6\nERFpScUmf7X8RUSaVlGTvOpt3hwe7L55c1jzR0SkkmV+kle9Xr1gwABYujR2JCIi6VSRyR9U9xcR\naY6Sv4hIBlV08lenr4hI4yqywxdgx45Q+6+rg4MPTiAwEZGUUodvgc6dYfRomD8/diQiIulTsckf\nVPcXEWlKxSd/1f1FRPZVsTV/gHXrYNQo2LQprPkjIlKJVPNvYMAA2H9/WLUqdiQiIulS0ckfYMwY\ndfqKiDRU8cn/xBPhpZdiRyEiki4Vn/zHjFHyFxFpqKI7fCE81WvgQHj//TD2X0Sk0qjDtxE9e8Jh\nh2mFTxGRQhWf/EGlHxGRhhJJ/mY22cyWm9kKM7uuiWNyZrbAzJaY2dNJXLe1NOJHRGRvRSd/M6sC\nbgcmASOBaWY2osExPYFZwNnufgxwfrHXbQuN+BER2VsSLf8xwEp3r3X37cBsYEqDY74CzHH3tQDu\n/k4C1221446DFSvg449LeVURkfRKIvkPBOoK9tfkPyt0JNDbzJ42s/lm9tUErttqXbvCMcfAggWl\nvKqISHqVavBjZ2A0MBHoDswzs3nu/kZjB8+cOXP3+1wuRy6XKzqA+tLPuHFFn0pEJKrq6mqqq6uL\nOkfR4/zNbCww090n5/evB9zdbyo45jqgm7vfkN+/G3jM3ec0cr5Ex/nXu+8+ePxxeOCBxE8tIhJV\nrHH+84HhZjbEzLoAU4FHGhzzMDDezDqZ2QHASUBNAtduNY34ERHZo+iyj7vvNLPpwFzCzeQed68x\ns8vDj/1Od19uZk8Ai4CdwJ3uvqzYa7fFUUeFpZ3ffRf69CnllUVE0qfil3coNHEifOc7MHlyh5xe\nRCQKLe/QApV+RESCTCV/TfYSEQkylfzr1/hJWaVLRKTkMpX8Bw2Cqiqoq2v5WBGRSpap5G+m0o+I\nCGQs+YOWdxYRgYwmf434EZGsy9Q4f4D33oPDD4fNm6FTpw67jIhIyWicfyv07g39+sHy5bEjERGJ\nJ3PJH1T6ERHJZPLXiB8RybpMJn+N+BGRrMtchy/AJ59A375hhc9u3Tr0UiIiHU4dvq20//5w5JGw\ncGHsSERE4shk8geVfkQk2zKd/DXiR0SyKrPJXyN+RCTLMtnhC7BjBxx8MKxZE15FRMqVOnzboHNn\nGD0aXnkldiQiIqWX2eQPKv2ISHZlOvlrxI+IZFXmk79G/IhIFmU6+R9+OGzdCmvXxo5ERKS0Mp38\nzdT6F5FsynTyByV/EcmmzCd/jfgRkSzK7CSvehs3hkXe3nsPqjJ/KxSRcqRJXu1wyCHQqxe88Ubs\nSERESieR5G9mk81suZmtMLPrmjnuRDPbbmbnJXHdpKj0IyJZU3TyN7Mq4HZgEjASmGZmI5o47kbg\niWKvmTRN9hKRrEmi5T8GWOnute6+HZgNTGnkuCuB3wAbE7hmojTiR0SyJonkPxCoK9hfk/9sNzMb\nAHzZ3X8KtKlTohRGj4ZFi+DTT2NHIiJSGp1LdJ1bgcK+gGZvADNnztz9PpfLkcvlOiSoej16wLBh\nsHgxfOELHXopEZGiVVdXU11dXdQ5ih7qaWZjgZnuPjm/fz3g7n5TwTGr6t8CfYGPgMvc/ZFGzlfS\noZ71/umfQsfvN75R8kuLiBQl1lDP+cBwMxtiZl2AqcBeSd3dh+W3oYS6/zcbS/wxacSPiGRJ0cnf\n3XcC04G5wFJgtrvXmNnlZnZZY7+k2Gt2BI34EZEsyfwM33rbt4fHOb79Nhx4YMkvLyLSbprhW4T9\n9oNRo+DVV2NHIiLS8ZT8C6j0IyJZoeRfQJO9RCQrlPwLaMSPiGSFkn+B4cPhgw/CMs8iIpVMyb9A\nVRWccIJKPyJS+ZT8GxgzBl58MXYUIiIdS8m/gZNPhnnzYkchItKxNMmrgXffhaFDw2MdO5dq2TsR\nkSJoklcC+vSBgQNhyZLYkYiIdBwl/0aMGwd/+UvsKEREOo6SfyNOOUXJX0Qqm5J/I9TyF5FKp+Tf\niCOPhI8/hjVrYkciItIxlPwbYRZKP88/HzsSEZGOoeTfBNX9RaSSKfk3Ydw4tfxFpHJpklcTtm4N\nY/43boTu3WNHIyLSNE3ySlC3buHJXlriWUQqkZJ/MzTkU0QqlZJ/M1T3F5FKpZp/MzZsgBEjwmJv\nVbpNikhKqeafsH79QqfvsmWxIxERSZaSfwtU+hGRSqTk3wJ1+opIJVLyb8H48fDss7GjEBFJlpJ/\nC44+GrZsgbq62JGIiCRHyb8FZjBhAjzzTOxIRESSk0jyN7PJZrbczFaY2XWN/PwrZrYwvz1nZscm\ncd1SUfIXkUpTdPI3syrgdmASMBKYZmYjGhy2Cpjg7scB/wbcVex1S+m005T8RaSyJNHyHwOsdPda\nd98OzAamFB7g7i+4+wf53ReAgQlct2SOPRbeeQfWr48diYhIMpJI/gOBwu7QNTSf3C8FHkvguiVT\nVQWnnqrWv4hUjs6lvJiZnQ5cAoxv7riZM2fufp/L5cjlch0aV2vUl36mTo0diYhkXXV1NdXV1UWd\no+i1fcxsLDDT3Sfn968H3N1vanDcKGAOMNnd32zmfKlZ26fQq6/ChRdCTU3sSERE9hZrbZ/5wHAz\nG2JmXYCpwCMNAhtMSPxfbS7xp9lxx4WF3lT3F5FKUHTyd/edwHRgLrAUmO3uNWZ2uZldlj/sB0Bv\n4A4zW2BmZfeIlE6dIJeDp56KHYmISPG0pHMbzJoVyj/33BM7EhGRPbSkcwebOBH+9CdI6b1JRKTV\nlPzbYMQI+PRTeOut2JGIiBRHyb8NzPa0/kVEypmSfxudcYaSv4iUP3X4tlFtLZx4Irz9tp7rKyLp\noA7fEhgyBHr1goULY0ciItJ+Sv7tMGkSPPFE7ChERNpPyb8dlPxFpNyp5t8OH30Ehx4alnro0SN2\nNCKSdar5l0j37jBmDDz9dOxIRETaR8m/nVT6EZFypuTfTkr+IlLOlPzbadSoUPtfsSJ2JCIibafk\n305mcPbZ8Ic/xI5ERKTtlPyLcPbZ8Pvfx45CRKTtNNSzCB99BP37Q10d9OwZOxoRySoN9Syx7t3h\n1FPV8Ssi5UfJv0gq/YhIOVLZp0h1dTB6dFjls1On2NGISBap7BPBYYeF7bnnYkciItJ6Sv4JOO88\neOih2FGIiLSeyj4JWLYszPitrdUDXkSk9FT2ieRznwure778cuxIRERaR8k/IX//9zBnTuwoRERa\nR8k/IfXJv8wqViKSUUr+CTn+eNi1S8/2FZHyoOSfEDP4x3+EBx6IHYmISMs02idBS5bAWWfBX/+q\nUT+l8uGH8NZbYVu1Kky2e+892Lw5vG7ZAjt2wM6d4dUdunXbs+2/P/TqBX36hK1v3z3v+/cPczgO\nOij271Kkee0Z7dM5oQtPBm4lfJO4x91vauSY24AzgY+Ai939tSSunSbHHAMHHxwmfE2YEDuayrNp\nE8yfH0ZV1W/vvw9Dh8KwYeG1f//wvnfvsPXoAZ07h61Tp/ANbds22LoVPvkEPv443CjefTdsixeH\n13feCc9orqsLv27QoHAjqH+t34YMgcGDw41EpJwU3fI3sypgBXAGsA6YD0x19+UFx5wJTHf3vzOz\nk4Afu/vYJs5Xti1/gBtvDC3/n/0sdiTlb/t2mDcPHn8cHnsstO5PPBFOOCFsX/hCSL7WpvZO27iH\nG8yaNeFGUP9av9XWhtdevUIsQ4bA4YfveV+/6duDdKT2tPyTSP5jgRnufmZ+/3rAC1v/ZvYz4Gl3\n/3V+vwbIufuGRs5X1sm/tjYkpXXroEuX2NGUn1274Nln4Re/gN/+NrTmzzwTJk+GsWNDCz5tdu0K\n3xJqa/fe/vrXPe+7dGn+5tC3b8fexKSyxSr7DATqCvbXAGNaOGZt/rN9kn+5GzIEjj4aHn0Uvvzl\n2NGUj9pauOsu+OUvQyv5oovgX/8VBg6MHVnLqqpCnAMHwimn7Ptz91BKanhT+POf93y2bVsoHw0a\nBP36Nb4dckjY9tuv5L/FzHHfe9h2Y+3Rhp+1tA/pahCmsB1V/i65BO69V8m/NV55BW65JTwT4aKL\n4OGH4bjjYkeVLLPQsu/bN3wrbMzf/garV4ey0saNsGFD2BYv3nv/nXega9fQt1S49ey553337qEj\nu6mta9fQj1FVFV6be79rV9h27mz9644d6dq2b2/7r3Hf95tYY9/MWjqmcP+II6Cmpu3/fjpKEsl/\nLTC4YH9Q/rOGxxzWwjG7zZw5c/f7XC5HLpcrNsaSuuAC+Pa3Q+lnwIDY0aTTn/4E//7v8MYbcPXV\noY8ky3Xxgw4KAwaOOab549zDCKb334cPPgivDbcPPww3jE8+aXzbtm1Poi5M2g3f79y550bQmtf6\n9/vtt6eTvZjtgAOSOU/h1trY0j5ar7q6murq6qLOkUTNvxPwOqHDdz3wEjDN3WsKjjkLuCLf4TsW\nuLVSO3zr/fM/h1En3/1u7EjS5ZVX4PrrQ6ljxoxwo1QZQ6Q4URZ2c/edwHRgLrAUmO3uNWZ2uZld\nlj/mUeAtM3sD+H/AN4u9btpdeincc09oSQm8+SZMnQrnnBOWwli6FC68UIlfJBZN8uog7jBqFNx2\nG5x+euxo4tm2DW66Kfw5XH01XHNNqEmLSHK0pHOKmME3vgG33x47knieeircABcsgFdfhe9/X4lf\nJC3U8u9AW7aEoZ+vvBLGdmfF5s1w1VXwzDPwk5/AuefGjkiksqnlnzI9esDFF8OsWbEjKZ25c0Nr\nv2fPUNdX4hdJJ7X8O1j9kgS1tZVd8vjoI7juujBO/9574Utfih2RSHao5Z9CQ4fCqafCfffFjqTj\nLFwIo0eHceeLFinxi5QDtfxLYN48mDYNVq6svKGN990H//IvcOutYeimiJRelIXdklaJyR9Ca3ja\nNPj612NHkoytW/d06s6ZAyNHxo5IJLtU9kmxH/wAfvjDsG5IuXvrLRg3LjwsZf58JX6RcqTkXyIT\nJoRVH8v9MY+PPhqWVv7qV+HBB+HAA2NHJCLtobJPCT3zTBj6WVNTfk9+2rkTbrghjOSZPRvGj48d\nkYjUU9kn5U47LSxXfNttsSNpm02bwgNVnn02TFhT4hcpf0r+JXbzzWHbtCl2JK3zwgthDfrRo+HJ\nJ8NDRUSk/KnsE8G3vhUeMPHTn8aOpGnucMcdodRz110wZUrsiESkKRrqWSY2bw4P7fj1r9NZQtmy\nBS67DJYtC8M4jzgidkQi0hzV/MtEr15hwbNLLw3j5dNk+XI46aTwqL9585T4RSqVkn8k550XWv83\n3BA7kj0efDAsRXHNNWFUz/77x45IRDqKyj4RbdgQOlLvvRcmTYoXxyefwLXXhhU5H3yw6YeMi0g6\nqexTZvr1g/vvh699DVavjhPD8uVh0tbmzeGBK0r8Itmg5B/ZaaeFVvc//EPoaC0Vd/jP/wxlnunT\nw8zjnj1Ld30RiUtlnxRwD52/dXXwP/8TOls70po14RGTq1fDr34Fxx7bsdcTkY6lsk+ZMoM774SD\nDgorf27b1jHXcQ9j9j//eRgzBl5+WYlfJKvU8k+RbdvCmvgbN8Jvfwt9+iR37mXL4Mor4cMPQwfz\nMcckd24RiUst/zLXtWsYbXPyyWGs/bx5xZ9z3Tq44grI5eCcc+D555X4RUTJP3WqquCmm+BHPwpz\nAa66Ct59t+3nqamBb34zJPpu3cL+1VdD587Jxywi5UfJP6XOPx+WLAkzgIcPDyNy5s0LSys35Y03\nwszh8eNh4kTo3TsM5bzllmRLSCJS/lTzLwPr14dF1h55JIzUGTkSBg8OLfqtW2Ht2lDT79QJJk8O\ni7CddVblPS9YRBqnhd0yYP16eP11qK2FTz8NN4BDDw03hP79w8ghEckWJX8RkQwq+WgfM+tlZnPN\n7HUze8LM9pkjamaDzOwpM1tqZovN7FvFXFNERIpXbIfv9cAf3f0o4Cngu40cswO41t1HAicDV5jZ\niCKvG1V1dXXsEFpFcSZLcSZLccZVbPKfAtyXf38f8OWGB7j72+7+Wv79FqAGGFjkdaMql38MijNZ\nijNZijOuYpP/Ie6+AUKSBw5p7mAzOxw4HnixyOuKiEgRWpzyY2ZPAoWP7TbAge83cniTPbVm1gP4\nDXBV/huAiIhEUtRoHzOrAXLuvsHMDgWedvejGzmuM/B74DF3/3EL59RQHxGRNmrraJ9iJ/s/AlwM\n3AR8DXi4iePuBZa1lPih7b8BERFpu2Jb/r2BB4HDgFrgAnd/38z6A3e5+9lmNg74M7CYUBZy4Hvu\n/njR0YuISLukbpKXiIh0vNQs7GZmk81suZmtMLPrYsfTmHKbsGZmVWb2qpk9EjuWpphZTzP7bzOr\nyf+5nhQ7pobM7BozW2Jmi8zsV2bWJXZM9czsHjPbYGaLCj5rcfJlCmK8Of93/pqZzTGzg2LGmI9p\nnzgLfvZtM9uVr3ZE1VScZnZl/s90sZnd2NJ5UpH8zawKuB2YBIwEpqV0Ili5TVi7ClgWO4gW/Bh4\nND9Q4DjCPJDUMLMBwJXAaHcfRegnmxo3qr38nPD/plBrJl+WUmMxzgVGuvvxwErixwiNx4mZDQK+\nRChtp8E+cZpZDjgHONbdjwX+o6WTpCL5A2OAle5e6+7bgdmECWSpUk4T1vL/YM8C7o4dS1Pyrb1T\n3f3nAO6+w93/FjmsxnQCuudHrR0ArIscz27u/hywucHHLU6+LKXGYnT3P7r7rvzuC8CgkgfWQBN/\nlgD/F/jfJQ6nSU3E+b+AG919R/6Yd1o6T1qS/0CgrmB/DSlNqvXKYMJa/T/YNHfqDAXeMbOf58tT\nd5rZ/rGDKuTu64BbgNXAWuB9d/9j3Kha1KbJlynwdeCx2EE0xszOBercfXHsWFpwJDDBzF4ws6fN\n7ISWfkFakn9ZSfuENTP7O2BD/luK5bc06gyMBma5+2jgY0LJIjXM7GBCS3oIMADoYWZfiRtVm6W2\nAWBm/wfY7u73x46loXxD5HvAjMKPI4XTks5AL3cfC3yHMAqzWWlJ/muBwQX7g/KfpU7+q/9vgP9y\n96bmNcQ2DjjXzFYBDwCnm9kvIsfUmDWEVtXL+f3fEG4GafJFYJW7v+fuO4GHgFMix9SSDWbWDyA/\n+XJj5HgaZWYXE0qTab2ZHgEcDiw0s7cIeekVM0vjN6k6wr9N3H0+sMvMmn1+X1qS/3xguJkNyY+k\nmEqYQJZGrZ6wFou7f8/dB7v7MMKf5VPuflHsuBrKlybqzOzI/EdnkL4O6tXAWDPrZmZGiDFVndLs\n++2ufvIlND/5spT2itHMJhPKkue6+7ZoUe1rd5zuvsTdD3X3Ye4+lNBY+by7p+Fm2vDv/HfARID8\n/6f93L3Zp3+nIvnnW1TTCSMAlgKz3T1t/8HIT1i7EJhoZgvyderJseMqc98CfmVmrxFG+/wwcjx7\ncfeXCN9IFgALCf/h7owaVAEzux94HjjSzFab2SXAjcCXzOx1ws2qxWF/EWL8CdADeDL//+iOmDFC\nk3EWclJQ9mkiznuBYWa2GLgfaLGxp0leIiIZlIqWv4iIlJaSv4hIBin5i4hkkJK/iEgGKfmLiGSQ\nkr+ISAYp+YuIZJCSv4hIBv1/sts9OwPgMpUAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.plot(ts, position, 'b');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's velocity as a function of time." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUVNW59/HvA+2AIAjIPIhKRBzBAZwSC5VEieMdctXk\nKr53ubJMTLyJ8Y0acyU317s0KzExmlfjEIO5iRpNFNSbMCwoEQVEAUVsEA0qY4PIEJGp6ef9Y1dL\n01TRw6muc07V77NWra4udp/zMP1q1z5772PujoiIVJZ2cRcgIiKlp/AXEalACn8RkQqk8BcRqUAK\nfxGRCqTwFxGpQEUJfzM738wWm9k7Zvb9PL8+xMxeMbNtZvbdYpxTRERaz6LO8zezdsA7wLnAKmAu\ncLm7L27Q5lDgMOBSYIO73x3ppCIiEkkxev4jgKXu/oG77wSeAC5p2MDdP3L314HaIpxPREQiKkb4\n9wOWN/h+Re41ERFJKF3wFRGpQFVFOMZKYGCD7/vnXmsVM9NmQyIiLeTu1pL2xej5zwUGm9lhZrY/\ncDkwcR/tmyzQ3RP9uP3222OvQXWqTtWpOusfrRG55+/uu8zsemAy4c3kEXevNrOvh1/2B82sF/Aa\ncDBQZ2Y3AMe4+ydRzy9tYPNmqK6GDz4Iz/fbD7p2hcGD4cgj4YAD4q5QRCIqxrAP7v5XYEij137d\n4HkNMKAY55I2snEjPPIIPPssLFgAQ4bAoEHQpQvU1sL69bB0KaxcCSNHwhe/CJdfDocdFnflItIK\nRQn/SpPJZOIuoVmaVefmzfDjH4fgHzMGbrsNvvAF6NAhf/tNm+Cll+D55+Hkk+G44+CGG+Dii6F9\n+7arMwFUZ3GpznhFXuRVbGbmSaupbE2YANdfD6NHhzeAfi2cobtjR/ikcPfd8NFH8B//AV/9aqvf\nBESkdcwMb+EFX4V/Jdq1K/Twn3gCfvtbOPvs6Md88cVwzPXr4Ze/hPPOi35MEWkWhb80bds2+MpX\nwnDP00/DoYcW79ju8NxzYRho5MjwiaBv3+IdX0Tyak34a5FXJdm+Hf7hH+Cgg2Dy5OIGP4BZGPtf\ntCjMDDrhBPjFL8IFYxFJFPX8K0VtLfzjP4Zpm088AVUluNa/eDF84xvh08Zjj4U3BBEpOvX8pbCb\nbgoh/PjjpQl+gKOPhqlTw5TQ00+HBx4IQ0MiEjv1/CvBY4/Bf/4nzJ0bFmvFoboa/vVfoUcP+M1v\noE+feOoQKUPq+cve3nwTbrwxTOuMK/gBhg6FWbNgxAgYPjysExCR2KjnX8527Ahhe8MNcM01cVez\n28yZ8LWvhYvDP/kJHHhg3BWJpJp6/rKnH/0IBg6EsWPjrmRPZ50F8+fD6tVhSmh1ddwViVQchX+5\nev31sGXDgw+GKZhJ07Ur/PGP8K1vhe0kHnxQF4NFSkjDPuWorg7OOAO+/vVkDfcUUl0NV1wRdgx9\n6CHo1i3uikRSRcM+Ejz2WPh69dXx1tFcQ4fC7NnQv3+4GPzSS3FXJFL21PMvN5s2hTCdMAFOPTXu\nalruhRfg3/4tfGr54Q9LtyZBJMW0t4+EzdVWroRHH427ktZbvRquugq2boXf/173DBBpgoZ9Kt26\ndXD//TBuXNyVRNOnD0yaFKaCnnoqPPVU3BWJlB31/MvJ974XtnC47764KymeV1+FK6+ETAbuuQc6\ndoy7IpHEUc+/kq1aFYZ6br017kqKa8SIsCZgx45w57DZs+OuSKQsqOdfLr7zHWjXDn72s7graTtP\nPhlWK19xBfzXf+lTgEiOev6VasMGGD8+vAGUs3/5F3jrrXDLyOOPDzuGikirKPzLwQMPwEUXhXny\n5e7QQ+F3v4Nf/SpMCf3qV2HFirirEkkdhX/abdsW7pn7ve/FXUlpXXABvP02HHEEDBsGd9wR/ixE\npFmKEv5mdr6ZLTazd8zs+wXa/NLMlprZAjMbVozzCmEe/PDhYRik0nTsCD/+cbhPweuvwzHHhLuU\n1dXFXZlI4kUOfzNrB9wHfAk4FrjCzI5u1OYC4Eh3/xzwdeCBqOcVwkZo990XLoJWssMPhz//GR5+\nONwz+MQTw83p9SYgUlAxev4jgKXu/oG77wSeAC5p1OYS4DEAd58DdDGzXkU4d2WbPRv+/ncYPTru\nSpLhnHPCDWPuugvuvDN8InrqKd1AXiSPYmyc0g9Y3uD7FYQ3hH21WZl7raYI569c998P110XpnhK\nYAZjxoRrAs89F24Wc+ON4c/p2mvDBWNpntpa2Lw57BfV+LFlS7jGsnVreNQ/r/+6fTvs2hU+fe3a\ntefzxq+ZhX/D7duHr/t6Xv99VVXrHlF+tjXHTfD/zUTumjWuwfYEmUyGTCYTWy2J9dFHIdx+/vO4\nK0kms7A9xMUXw7x5cO+98LnPwSWXhBlCo0ZV3qZxtbVhMeCqVbB2bXisW7f313XrYOPGEOIHHwxd\nuuz96NQp3IGtQ4fw6N599/MDD4QDDgiBWB/cjZ83/Aq73wgavkE0ft7wTWPXrvD7aeqxbVv+15v7\n8y19NDzuzp3h32H9G8FRR8GCBUX5q8xms2Sz2UjHiLzIy8xOA8a5+/m5728G3N3vatDmAWC6uz+Z\n+34xcLa779Xz1yKvZvrpT8Oc99/+Nu5K0mPdujBN9PHHYfly+Od/DmsHTjst/W8EdXXh9/fhh+H3\n1vBR/9ratdCjB/TrBz17hkePHnt+7dkzfDo65JAQ8AnuuaZCXd3uN4O6uvBn2gZi2dXTzNoDS4Bz\ngdXAq8AV7l7doM0Y4Jvu/uXcm8Uv3P20AsdT+DfFPcxsefhhOPPMuKtJp3ffDTODnnoqhOOoUXDe\neeH6yeDBybr7mXvoiRcK9eXLw06uBx8MAwbs+Rg4cPfzvn1hv/3i/t1IG4htS2czOx+4h3AB+RF3\nv9PMvk74BPBgrs19wPnAFuAad59X4FgK/6bMmRNugP7OO8kKqbRasyasFq5/bN0aLhafdFL4OnQo\nDBoUesPFVlsLNTVhKGblyvBo+HzlyhDusGeQNw72/v3hoIOKX5+kgvbzrxTf+Eb46P6DH8RdSXla\nvTpsJjd/frhesHQpLFsWhoYGDQp/9l277n507rz3Rb7t23c/tm0LPfcNG8Lj44/3/Nq9e+iV9+sX\nHvXP678OGBDG2fVGLwUo/CvBtm0hEObPDz0/KQ33ENbLloU3h/og37AhzIhpeKFv165wwbPh45BD\nwr2Ju3bd/bVr1zC+rqEYiag14Z/yq1wVaOLEMBSh4C8ts9BD79497kpEikKX8tNm/Phwi0MRkQg0\n7JMm69eHjcxWrAgzO0RE0H7+5e+ZZ8JURAW/iESk8E+TJ5+Eyy+PuwoRKQMa9kmLmhoYMiTMAdd8\nbhFpQMM+5exPf4Ivf1nBLyJFkczwV89/b08+GfahEREpgmQO+3zySbhLkwRr14YdAdesCTsmiog0\nUD7DPh9/HHcFyTJxInzpSwp+ESkahX8aPPMMXHpp3FWISBlJZvhv2BB3BcmxeTO89FK4O5WISJEk\nM/zV89/tL38Je/Z36RJ3JSJSRhT+SffMM3DZZXFXISJlJpnhr2GfYMcOmDQp3IdWRKSIkhn+6vkH\nM2eGm4737h13JSJSZhT+SfbCC3DhhXFXISJlKJnhr2Gf4IUXwpYOIiJFlszwV88f3n0XNm0Kd+0S\nESkyhX9S1ff62yXzr0hE0i2ZyaLw15CPiLSpZG7s1rlzGPKoVFu2hBk+q1bprl0i0qSSb+xmZl3N\nbLKZLTGzSWaWdxmqmT1iZjVm9mazDvzJJ7BzZ5TS0m3GDDjpJAW/iLSZqMM+NwNT3X0IMA24pUC7\nR4EvNfuohxwCGzdGLC3Fpk4N9+oVEWkjUcP/EmB87vl4IO/Wk+4+E2j+/M1u3Sp73H/KFIW/iLSp\nqOHf091rANx9DdAzeklA9+6wfn1RDpU6a9bA8uVw8slxVyIiZayqqQZmNgXo1fAlwIHb8jQvytXj\ncR9/DPfeC5Mnk8lkyGQyxThsOkydCqNGQVWTfzUiUqGy2SzZbDbSMSLN9jGzaiDj7jVm1huY7u5D\nC7Q9DHjO3U9o4pjuV10FmQxcc02ra0utsWNh5Ei47rq4KxGRlIjjNo4TgbG551cDE/bR1nKPph16\nKHz0UaTCUsld4/0iUhJRw/8uYLSZLQHOBe4EMLM+ZvZ8fSMz+wPwCnCUmX1oZvvu0ldq+FdXw377\nwZFHxl2JiJS5SAPL7v4xcF6e11cDFzb4/soWHfjQQ+G996KUlk71UzytRZ/eRERaLJnbO3TvXpk9\nfw35iEiJJDP8K3HYZ+fOsLL3nHPirkREKkByw7/S5vnPng2DB4ffu4hIG0tu+Fdaz19bOohICSUz\n/Lt1C3fz2rUr7kpKR+P9IlJCyQz/qqqwo2WlbO62aRMsXAhnnhl3JSJSIZIZ/lBZ4/4zZ8KIEXDg\ngXFXIiIVItnhXynj/i++CGefHXcVIlJBFP5JMGOGwl9ESiq54V8pC70++QTeeisM+4iIlEhyw79S\nxvxnzYLhw6FDh7grEZEKkuzwr4Sev8b7RSQGCv+4zZgBX/hC3FWISIVJbvhXwq0ct26FefPgjDPi\nrkREKkxyw78Sev6vvgrHHgudOsVdiYhUGIV/nDTFU0RiovCP04svarxfRGIR6QbubcHM3N2htjZs\nd7B9O7RvH3dZxbdjR7iusXw5HHJI3NWISIrFcQP3tlNVBZ07l+/mbq+/HvbvV/CLSAySG/5Q3kM/\nGu8XkRgp/OOi8X4RiVGyw79HD1i3Lu4qim/XLnjlFfj85+OuREQqVLLDv2dPWLs27iqKb8EC6Ncv\nvLmJiMQgUvibWVczm2xmS8xskpl1ydOmv5lNM7NFZrbQzL7d7BP07Ak1NVFKTKaZM+Gss+KuQkQq\nWNSe/83AVHcfAkwDbsnTphb4rrsfC5wOfNPMjm7W0Xv1Ks+e/6xZ2tJBRGIVNfwvAcbnno8HLm3c\nwN3XuPuC3PNPgGqgX7OOXq7DPrNmwemnx12FiFSwqOHf091rIIQ80HNfjc1sEDAMmNO8o5fhsM/K\nlbBlC3zuc3FXIiIVrKqpBmY2BejV8CXAgdvyNC+4XNjMOgFPAzfkPgEUNG7cuPBk7Voyy5aRaarI\nNKnv9VuLFuOJiHwmm82SzWYjHSPS9g5mVg1k3L3GzHoD0919aJ52VcDzwF/c/Z4mjumf1bRuHQwd\nWl5z/W+8MWzrcOutcVciImUiju0dJgJjc8+vBiYUaPcb4O2mgn8v3brBpk2wc2erC0wcXewVkQSI\n2vPvBvwRGAB8AHzF3TeaWR/gIXe/0MzOBGYACwnDQg7c6u5/LXBM36Om3r1h/nzo06fVdSbG9u3h\nDW3tWujYMe5qRKRMtKbn3+SY/764+8fAeXleXw1cmHv+MtD6bTnrZ/yUQ/jPmwdDhij4RSR2yV7h\nC+U140dTPEUkIZIf/uW00Evj/SKSEMkP/3JZ6OUeNnNTz19EEiAd4V8Owz7Ll4e7kx1+eNyViIik\nIPzLZdhHi7tEJEGSH/7lMuzzyisa7xeRxEhH+JfDsI9m+ohIgiQ//Mth2GfrVli0CE45Je5KRESA\nNIR/jx6h5x9hJXLsXn8djjkGOnSIuxIRESAN4X/QQbD//rB5c9yVtJ7G+0UkYZIf/pD+oR+N94tI\nwqQj/NN80ddd4S8iiaPwb2vLlkH79jBwYNyViIh8Jh3h36dPesNfi7tEJIHSE/6rVsVdRetoyEdE\nEigd4d+3r8JfRKSI0hH+ffrA6tVxV9FyW7bA4sVw0klxVyIisod0hH9ae/5z58IJJ8CBB8ZdiYjI\nHtIT/mns+WvIR0QSKh3h36MHbNwIO3bEXUnL6OYtIpJQ6Qj/du3SN9ffHWbPVviLSCKlI/whfdM9\n3303bOTWv3/clYiI7CU94Z+2cX+N94tIgkUKfzPramaTzWyJmU0ysy552hxgZnPMbL6ZLTSz21t1\nsrT1/LWTp4gkWNSe/83AVHcfAkwDbmncwN23A6PcfTgwDLjAzEa0+Expm+6pnr+IJFjU8L8EGJ97\nPh64NF8jd/809/QAoApo+Z1Z0jTss3kzvPceDBsWdyUiInlFDf+e7l4D4O5rgJ75GplZOzObD6wB\nprj73BafKU3DPq++GoJ///3jrkREJK+qphqY2RSgV8OXCD332/I0z9ujd/c6YLiZdQaeNbNj3P3t\nQuccN27cZ88zmQyZTCZdPf9ZszTeLyJtJpvNks1mIx3DPMK9cc2sGsi4e42Z9Qamu/vQJn7mh8AW\nd7+7wK973prWrAlbJaThjl5jxsC118Jll8VdiYhUADPD3Vu0b3zUYZ+JwNjc86uBCXmKOrR+FpCZ\ndQBGA4tbfKb6Vb47d7a62JKoq9PiLhFJvKjhfxcw2syWAOcCdwKYWR8zez7Xpg8w3cwWAHOASe7+\nvy0+U/v2YZXvmjURS25jS5ZAly7Qu3fclYiIFNTkmP++uPvHwHl5Xl8NXJh7vhAozp7G9Rd9Bwwo\nyuHahMb7RSQF0rPCF9Jx0Vfz+0UkBdIV/mmY7qnwF5EUSFf4DxgAK1bEXUVhGzfC+++HWUkiIgmW\nvvD/8MO4qyhszhw45RTYb7+4KxER2ad0hf/AgckOfw35iEhKKPyLSXfuEpGUiLTCty0UXOELsH07\ndO4Mn34a5v0nSV0ddOsGS5eGBWkiIiUSxwrf0jrgAOjePZkLvd5+O4S+gl9EUiBd4Q/JHfrReL+I\npEj6wj+pM3403i8iKZK+8E9yz1/bOohISij8i2H9+rDy+Ljj4q5ERKRZFP7FMHs2nHpq8mYgiYgU\nkM7wX7487ir2pIu9IpIy6Qv/JF7w1Xi/iKRMuhZ5AbjDQQfBRx9Bx46lK6yQ2tqwuOv998NXEZES\nK/9FXgBmofeflKGfN94I9Sj4RSRF0hf+kKyLvjNnwllnxV2FiEiLKPyjevllhb+IpE46wz8pF33d\n1fMXkVRKZ/gfcQT87W9xVwHLloVrEIMGxV2JiEiLpDP8jzwS3nsv7ipCr//MM8MbgIhIiij8o9CQ\nj4ikVKTwN7OuZjbZzJaY2SQz67KPtu3MbJ6ZTYxyTgB694YtW2Dz5siHikQXe0UkpaL2/G8Gprr7\nEGAacMs+2t4AvB3xfIFZ/L3/9evDWoMTToivBhGRVooa/pcA43PPxwOX5mtkZv2BMcDDEc+32+DB\n8O67RTtci73yCowcCVVV8dUgItJKUcO/p7vXALj7GqBngXY/B24CireXRNw9f433i0iKNdltNbMp\nQK+GLxFC/LY8zfcKdzP7MlDj7gvMLJP7+X0aN27cZ88zmQyZTGbvRoMHw9y5TR2q7bz8MvzoR/Gd\nX0QqVjabJZvNRjpGpI3dzKwayLh7jZn1Bqa7+9BGbf4b+BpQC3QADgb+7O5XFTjmvjd2qzd1Ktxx\nB0yf3ur6W23btnAj+Zoa6NSp9OcXEWkgjo3dJgJjc8+vBiY0buDut7r7QHc/ArgcmFYo+FskzmGf\n116DY45R8ItIakUN/7uA0Wa2BDgXuBPAzPqY2fNRi9unAQNg7drQCy+1+sVdIiIplb79/Bs66iiY\nMAGGDm26bTFddBFcfTX80z+V9rwiInlUxn7+DcUx3bOuLlzsVc9fRFIs3eEfx7j/okXhxi19+pT2\nvCIiRZTu8I+j5z99OowaVdpziogUWbrDf8gQWLy4tOfMZhX+IpJ66b7gu3w5jBgBq1e3bVH16uqg\nRw94803o16805xQRaULlXfDt3x8+/TRsslYKCxeGxV0KfhFJuXSHv1lYbLVoUWnOp/F+ESkT6Q5/\ngOOOK134Z7OQb58hEZGUSX/4H3tsacJ/1y6YMUPhLyJloTzC/6232v48b7wBPXtqfr+IlIXyCP9S\n9PynTIEvfrHtzyMiUgLpD/8+faC2Nmzy1pYmT4bRo9v2HCIiJZL+8Ddr+4u+W7bAnDka7xeRspH+\n8Ie2H/efMQNOPhkOPrjtziEiUkLlE/5t2fPXeL+IlJnyCP/jjw+zcdqKxvtFpMyke2+fen//O/Tu\nDRs2wP77F7eglSvhhBPCBeX27Yt7bBGRIqi8vX3qHXwwHHFE2Hun2P76VzjvPAW/iJSV8gh/CLt7\nzplT/OM+91y4baOISBkpn/AfObL44b91K0ybBmPGFPe4IiIxK5/wHzECXn21uMecNg2GDw+3bRQR\nKSPlE/7HHRdu7rJxY/GOOXGihnxEpCyVT/hXVYVe+muvFed47vD88wp/ESlLkcLfzLqa2WQzW2Jm\nk8ysS4F275vZG2Y238yKPDbTQDHH/efNg44dw32CRUTKTNSe/83AVHcfAkwDbinQrg7IuPtwdx8R\n8ZyFFTP8n34aLrusOMcSEUmYSIu8zGwxcLa715hZbyDr7kfnabcMOMXdm7zZbqsWedVbtSqM/a9d\nG4aBWssdDj8cnn0Whg1r/XFEREogjkVePd29BsDd1wA9C7RzYIqZzTWzayOes7C+fWHAgOizfmbN\ngg4d4MQTi1OXiEjCNNk9NrMpQK+GLxHC/LY8zQt12c9099Vm1oPwJlDt7jMLnXPcuHGfPc9kMmRa\nspXy+eeHVblnnNH8n2ns8cfhyivDdtEiIgmTzWbJZrORjhF12KeaMJZfP+wz3d2HNvEztwN/d/e7\nC/x664d9INxk/aabYO7c1v18bS307w8zZ8Lgwa2vQ0SkROIY9pkIjM09vxqYkKeog8ysU+55R+CL\nQNttvn/GGbB0aevv7DV9ehg6UvCLSBmLGv53AaPNbAlwLnAngJn1MbPnc216ATPNbD4wG3jO3SdH\nPG9h++8Po0aFbZhb46GH4KqriluTiEjClMeWzo09+GC4+9b//E/Lfm7FirB98/vvQ+fO0WoQESmR\nyt3SubELLggXfbdta9nP/frX4UKvgl9Eylx5hv+AAeGeu08/3fyf2b49DPlcf33b1SUikhDlGf4A\n110H99/f/PZPPhmGfI7ea42aiEjZKc8xfwhTNg8/HF54IYT6vmzfDsccE64VnHtu9HOLiJSQxvwb\nqqqCa69tXu//l78M4a/gF5EKUb49f9i918+bb4aFW/msXRuC/+WXtYOniKSSev6N9e0L3/0ujB0L\ndXV7/7o7/Pu/w9e+puAXkYpS3uEPcPPN8OmncO+9e77uHraBWLYM7rgjntpERGISYd/jlKiqgt/9\nDk47DTZsgO98B9avD28GU6aExWAdO8ZdpYhISZV/zx/gyCPDTV7efz+M/Z92Wvg0MHmybs4uIhWp\nvC/45rNuHXTtGu1mLyIiCdKaC76VF/4iImVGs31ERKRZFP4iIhVI4S8iUoEU/iIiFUjhLyJSgRT+\nIiIVSOEvIlKBFP4iIhVI4S8iUoEU/iIiFShS+JtZVzObbGZLzGySmXUp0K6LmT1lZtVmtsjMRkY5\nr4iIRBO1538zMNXdhwDTgFsKtLsH+F93HwqcCFRHPG+sstls3CU0i+osLtVZXKozXlHD/xJgfO75\neODSxg3MrDPweXd/FMDda919c8Tzxiot/xhUZ3GpzuJSnfGKGv493b0GwN3XAD3ztDkc+MjMHjWz\neWb2oJl1iHheERGJoMnwN7MpZvZmg8fC3NeL8zTPtxdzFXAS8Ct3Pwn4lDBcJCIiMYm0n7+ZVQMZ\nd68xs97A9Ny4fsM2vYBZ7n5E7vuzgO+7+0UFjqnN/EVEWqil+/lHvZ3VRGAscBdwNTAhT0E1Zrbc\nzI5y93eAc4G3Cx2wpb8BERFpuag9/27AH4EBwAfAV9x9o5n1AR5y9wtz7U4EHgb2A/4GXOPum6IW\nLyIirZO42ziKiEjbS8wKXzM738wWm9k7Zvb9uOvJx8z6m9m03EK1hWb27bhr2hcza5ebYTUx7loK\nScMCQDP7jpm9lZvo8Hsz2z/umuqZ2SNmVmNmbzZ4rVmLL2Ou8Se5v/MFZvan3JTwWOWrs8Gv3Whm\ndbnRjlgVqtPMvpX7M11oZnc2dZxEhL+ZtQPuA74EHAtcYWZHx1tVXrXAd939WOB04JsJrbPeDezj\n+kpCJHoBoJn1Bb4FnOTuJxCuk10eb1V7eJTw/6ah5i6+LJV8NU4GjnX3YcBS4q8R8teJmfUHRhOG\ntpNgrzrNLANcBBzv7scDP23qIIkIf2AEsNTdP3D3ncAThAVkieLua9x9Qe75J4Sg6hdvVfnl/sGO\nIVxrSaQULQBsD3Q0syrgIGBVzPV8xt1nAhsavdzk4stSyleju09197rct7OB/iUvrJECf5YAPwdu\nKnE5BRWo8zrgTnevzbX5qKnjJCX8+wHLG3y/goSGaj0zGwQMA+bEW0lB9f9gk3xRJ/ELAN19FfAz\n4ENgJbDR3afGW1WTmrP4Mkn+D/CXuIvIJ7eeabm7L4y7liYcBXzBzGab2XQzO6WpH0hK+KeKmXUC\nngZuyH0CSBQz+zJQk/uUYrlHEiV+AaCZHULoSR8G9AU6mdmV8VbVYontAJjZD4Cd7v6HuGtpLNcR\nuRW4veHLMZXTlCqgq7ufBvxfwizMfUpK+K8EBjb4vn/utcTJffR/Gvidu++1riEhzgQuNrO/AY8D\no8zssZhrymcFoVf1Wu77pwlvBklyHvA3d//Y3XcBfwbOiLmmptTkFleSW3y5NuZ68jKzsYShyaS+\nmR4JDALeMLNlhFx63cyS+ElqOeHfJu4+F6gzs+77+oGkhP9cYLCZHZabSXE5YQFZEv0GeNvd74m7\nkELc/VZ3H5hbVX05MM3dr4q7rsZyQxPLzeyo3Ev7XAAYkw+B08zsQDMzQo2JuijN3p/u6hdfQoHF\nlzHYo0YzO58wLHmxu2+Praq9fVanu7/l7r3d/Qh3P5zQWRnu7kl4M238d/4scA5A7v/Tfu6+fl8H\nSET453pU1xNmACwCnnD3pP0Hw8zOBL4KnGNm83Pj1OfHXVfKfRv4vZktIMz2+e+Y69mDu79K+EQy\nH3iD8B/uwViLasDM/gC8AhxlZh+a2TXAncBoM1tCeLNqctpfDDXeC3QCpuT+H/2/OGuEgnU25CRg\n2KdAnb8BjjCzhcAfgCY7e1rkJSJSgRLR8xcRkdJS+IuIVCCFv4hIBVL4i4hUIIW/iEgFUviLiFQg\nhb+ISAUpU21BAAAADElEQVRS+IuIVKD/DyWVce/wYnQSAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.plot(ts, velocity, 'r');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Same example, now with units\n", "\n", "I'll import Pint and instantiate a registry." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import pint\n", "UNITS = pint.UnitRegistry()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For reasons we'll see later, I'm going to monkey patch the Unit and Quantity classes so they don't cache dimensionality information." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true }, "outputs": [], "source": [ "@property\n", "def dimensionality(self):\n", " \"\"\"Unit's dimensionality (e.g. {length: 1, time: -1})\n", " \"\"\"\n", " dim = self._REGISTRY._get_dimensionality(self._units)\n", " return dim\n", "\n", "pint.unit._Unit.dimensionality = dimensionality\n", "pint.quantity._Quantity.dimensionality = dimensionality" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here are the units we'll use for this problem." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [], "source": [ "kg = UNITS.kilogram\n", "s = UNITS.second\n", "N = UNITS.newton\n", "m = UNITS.meter" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here are their dimensions." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(,\n", " ,\n", " ,\n", " )" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "UNITS.kilogram.dimensionality, UNITS.second.dimensionality, UNITS.newton.dimensionality, UNITS.meter.dimensionality" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now I'll define the parameters of the problem with units." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mass = 1 * (kg)\n", "k = 1 * (N / m) # spring constant\n", "c = 1 * (N / (m/s)) # dynamic friction\n", "\n", "# times\n", "ts = np.arange(0, 15.0, 0.1) * s" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here are their dimensions:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(,\n", " ,\n", " ,\n", " )" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mass.dimensionality, k.dimensionality, c.dimensionality, ts.dimensionality" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's a version of the slope function with the parameters." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def slope_func(Y, t):\n", " y1, y2 = Y\n", " y1p = y2\n", " y2p = -(c*y2 + k*y1) / mass\n", " return y1p, y2p" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And the initial conditions with units." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[, ]" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y_init = [1 * (m), 0 * (m/s)]\n", "y_init" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now when we call the slope function the results have the right units." ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "(, )" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "slope_func(y_init, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But when we call `odeint`, we get an error, presumably because it was not expecting things with units." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "ename": "ValueError", "evalue": "setting an array element with a sequence.", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0masol\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0modeint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mslope_func\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0my_init\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[1;32m/home/downey/anaconda2/lib/python2.7/site-packages/scipy/integrate/odepack.pyc\u001b[0m in \u001b[0;36modeint\u001b[1;34m(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg)\u001b[0m\n\u001b[0;32m 213\u001b[0m output = _odepack.odeint(func, y0, t, args, Dfun, col_deriv, ml, mu,\n\u001b[0;32m 214\u001b[0m \u001b[0mfull_output\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrtol\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0matol\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtcrit\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mh0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhmax\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhmin\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 215\u001b[1;33m ixpr, mxstep, mxhnil, mxordn, mxords)\n\u001b[0m\u001b[0;32m 216\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0moutput\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m<\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 217\u001b[0m \u001b[0mwarning_msg\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_msgs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0moutput\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;34m\" Run with full_output = 1 to get quantitative information.\"\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;31mValueError\u001b[0m: setting an array element with a sequence." ] } ], "source": [ "asol = odeint(slope_func, y_init, ts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Turning off the units temporarily\n", "\n", "We can make everything behave as if it is dimensionless by replacing `UNITS._get_dimensionality` with a function that always returns an empty dictionary." ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ ">" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "saved_get_dimensionality_method = UNITS._get_dimensionality\n", "UNITS._get_dimensionality = lambda self: {}\n", "saved_get_dimensionality_method" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now all units have no dimensionality." ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "({}, {}, {}, {})" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "UNITS.kilogram.dimensionality, UNITS.second.dimensionality, UNITS.newton.dimensionality, UNITS.meter.dimensionality" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And all quantities have no dimensionality." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "({}, {}, {}, {})" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mass.dimensionality, k.dimensionality, c.dimensionality, ts.dimensionality" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can still call the slope function, and it still works." ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "(, )" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "slope_func(y_init, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now we can call `odeint`, too!" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": true }, "outputs": [], "source": [ "asol = odeint(slope_func, y_init, ts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The results are the same as what we saw before." ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": true }, "outputs": [], "source": [ "position, velocity = asol.transpose()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's position as a function of time." ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG9ZJREFUeJzt3Xm0VOWZ7/Hvc0BAQZEhIoMgSBSDoiGKKIglJgu0VbLs\n1oa4YjRta27EOKRvNLnJAnt1J+pq1zVGTK5T2nSixA4m2okDJnqMRlRUZDwIijkcBgEVjaggw3P/\neOtAcThz7VPvrtq/z1p7Ve06m70fpme/9bzDNndHRESypSp2ACIiUnpK/iIiGaTkLyKSQUr+IiIZ\npOQvIpJBSv4iIhmUSPI3s3vMbIOZLWrmmNvMbKWZvWZmxydxXRERaZ+kWv4/ByY19UMzOxM4wt0/\nC1wO/Cyh64qISDskkvzd/TlgczOHTAF+kT/2RaCnmfVL4toiItJ2par5DwTqCvbX5j8TEZEI1OEr\nIpJBnUt0nbXAYQX7g/Kf7cPMtNiQiEgbubu15fgkW/6W3xrzCHARgJmNBd539w1NnWj+fKdvX2fR\nIsc9fduMGTOix6A4FafiVJz1W3sk0vI3s/uBHNDHzFYDM4AugLv7ne7+qJmdZWZvAB8BlzR3vhNO\ngFtvhSlT4OWXoXfvJKIUEZF6iSR/d/9KK46Z3pZzXnghPPMM3Hgj3Hxz+2MTEZF9pbrDd8YMuPtu\nWL8+diR7y+VysUNoFcWZLMWZLMUZl7W3XtRRzMwLY7r2WtixA267LWJQIiIpZmZ4Gzt8U5/8N26E\nESPgtddg8OCIgYmIpFR7kn+qyz4AhxwCF18Md9wROxIRkcqR+pY/QE0NTJwIq1fDfvtFCkxEJKUq\nsuUPcPTRMGwYPPpo7EhERCpDWSR/gEsvDSN/RESkeGVR9gHYsgUOOwyWLIGBWhJORGS3ii37APTo\nAeefD7/8ZexIRETKX9kkf4ALLoA5c2JHISJS/sqm7AOwfTsceigsWKAx/yIi9Sq67ANhmOe558Lv\nfhc7EhGR8lZWyR/gvPPgoYdiRyEiUt7KquwDsHVrKP2sWBFm/4qIZF3Fl30AunWDyZPh4YdjRyIi\nUr7KLvlDqPv/4Q+xoxARKV9lV/YB2LQJhg8Pr126lCgwEZGUykTZB+Azn4GjjoLnn48diYhIeSrL\n5A+h7v/447GjEBEpT2Wb/CdNgieeiB2FiEh5KsuaP4RHO37mM7BsGfTvX4LARERSKjM1f4DOneGL\nX4S5c2NHIiJSfso2+YNKPyIi7VW2ZR+AVatg/HhYuxasTV94REQqR6bKPgBDh4bF3lasiB2JiEh5\nKevkbwa5HDz9dOxIRETKS1knf4DTT4fq6thRiIiUl7Ku+QPU1sJJJ8H69ar7i0g2Rav5m9lkM1tu\nZivM7LpGfn6QmT1iZq+Z2WIzuziJ6wIMGQIHHAA1NUmdUUSk8hWd/M2sCrgdmASMBKaZ2YgGh10B\nLHX344HTgVvMrHOx166nur+ISNsk0fIfA6x091p33w7MBqY0OMaBA/PvDwTedfcdCVwbCHV/JX8R\nkdZLIvkPBOoK9tfkPyt0O/A5M1sHLASuSuC6u02YAM89BynrvhARSa1SjfaZBCxw9wHA54FZZtYj\nqZMPHhzG+7/5ZlJnFBGpbEnU3dcCgwv2B+U/K3QJ8CMAd3/TzN4CRgAvN3bCmTNn7n6fy+XI5XLN\nBmAWZvo+91x4yIuISCWrrq6musgx7kUP9TSzTsDrwBnAeuAlYJq71xQcMwvY6O43mFk/QtI/zt3f\na+R8bRrqWW/WLFiwAO6+u52/ERGRMhVlqKe77wSmA3OBpcBsd68xs8vN7LL8Yf8GnGJmi4Ange80\nlviLUd/yFxGRlpX9JK96O3dCnz6wcmVY519EJCsyt7BboU6d4OST4S9/iR2JiEj6VUzyB5V+RERa\nS8lfRCSDKqbmD/Dxx6He/9570LVrwoGJiKRUpmv+EBZ4O+qoMORTRESaVlHJH2DsWHjhhdhRiIik\nW0Um/xdfjB2FiEi6VWTyV8tfRKR5FZf8P/tZ+OADePvt2JGIiKRXxSV/s/BYR5V+RESaVnHJH1T6\nERFpScUmf7X8RUSaVlGTvOpt3hwe7L55c1jzR0SkkmV+kle9Xr1gwABYujR2JCIi6VSRyR9U9xcR\naY6Sv4hIBlV08lenr4hI4yqywxdgx45Q+6+rg4MPTiAwEZGUUodvgc6dYfRomD8/diQiIulTsckf\nVPcXEWlKxSd/1f1FRPZVsTV/gHXrYNQo2LQprPkjIlKJVPNvYMAA2H9/WLUqdiQiIulS0ckfYMwY\ndfqKiDRU8cn/xBPhpZdiRyEiki4Vn/zHjFHyFxFpqKI7fCE81WvgQHj//TD2X0Sk0qjDtxE9e8Jh\nh2mFTxGRQhWf/EGlHxGRhhJJ/mY22cyWm9kKM7uuiWNyZrbAzJaY2dNJXLe1NOJHRGRvRSd/M6sC\nbgcmASOBaWY2osExPYFZwNnufgxwfrHXbQuN+BER2VsSLf8xwEp3r3X37cBsYEqDY74CzHH3tQDu\n/k4C1221446DFSvg449LeVURkfRKIvkPBOoK9tfkPyt0JNDbzJ42s/lm9tUErttqXbvCMcfAggWl\nvKqISHqVavBjZ2A0MBHoDswzs3nu/kZjB8+cOXP3+1wuRy6XKzqA+tLPuHFFn0pEJKrq6mqqq6uL\nOkfR4/zNbCww090n5/evB9zdbyo45jqgm7vfkN+/G3jM3ec0cr5Ex/nXu+8+ePxxeOCBxE8tIhJV\nrHH+84HhZjbEzLoAU4FHGhzzMDDezDqZ2QHASUBNAtduNY34ERHZo+iyj7vvNLPpwFzCzeQed68x\ns8vDj/1Od19uZk8Ai4CdwJ3uvqzYa7fFUUeFpZ3ffRf69CnllUVE0qfil3coNHEifOc7MHlyh5xe\nRCQKLe/QApV+RESCTCV/TfYSEQkylfzr1/hJWaVLRKTkMpX8Bw2Cqiqoq2v5WBGRSpap5G+m0o+I\nCGQs+YOWdxYRgYwmf434EZGsy9Q4f4D33oPDD4fNm6FTpw67jIhIyWicfyv07g39+sHy5bEjERGJ\nJ3PJH1T6ERHJZPLXiB8RybpMJn+N+BGRrMtchy/AJ59A375hhc9u3Tr0UiIiHU4dvq20//5w5JGw\ncGHsSERE4shk8geVfkQk2zKd/DXiR0SyKrPJXyN+RCTLMtnhC7BjBxx8MKxZE15FRMqVOnzboHNn\nGD0aXnkldiQiIqWX2eQPKv2ISHZlOvlrxI+IZFXmk79G/IhIFmU6+R9+OGzdCmvXxo5ERKS0Mp38\nzdT6F5FsynTyByV/EcmmzCd/jfgRkSzK7CSvehs3hkXe3nsPqjJ/KxSRcqRJXu1wyCHQqxe88Ubs\nSERESieR5G9mk81suZmtMLPrmjnuRDPbbmbnJXHdpKj0IyJZU3TyN7Mq4HZgEjASmGZmI5o47kbg\niWKvmTRN9hKRrEmi5T8GWOnute6+HZgNTGnkuCuB3wAbE7hmojTiR0SyJonkPxCoK9hfk/9sNzMb\nAHzZ3X8KtKlTohRGj4ZFi+DTT2NHIiJSGp1LdJ1bgcK+gGZvADNnztz9PpfLkcvlOiSoej16wLBh\nsHgxfOELHXopEZGiVVdXU11dXdQ5ih7qaWZjgZnuPjm/fz3g7n5TwTGr6t8CfYGPgMvc/ZFGzlfS\noZ71/umfQsfvN75R8kuLiBQl1lDP+cBwMxtiZl2AqcBeSd3dh+W3oYS6/zcbS/wxacSPiGRJ0cnf\n3XcC04G5wFJgtrvXmNnlZnZZY7+k2Gt2BI34EZEsyfwM33rbt4fHOb79Nhx4YMkvLyLSbprhW4T9\n9oNRo+DVV2NHIiLS8ZT8C6j0IyJZoeRfQJO9RCQrlPwLaMSPiGSFkn+B4cPhgw/CMs8iIpVMyb9A\nVRWccIJKPyJS+ZT8GxgzBl58MXYUIiIdS8m/gZNPhnnzYkchItKxNMmrgXffhaFDw2MdO5dq2TsR\nkSJoklcC+vSBgQNhyZLYkYiIdBwl/0aMGwd/+UvsKEREOo6SfyNOOUXJX0Qqm5J/I9TyF5FKp+Tf\niCOPhI8/hjVrYkciItIxlPwbYRZKP88/HzsSEZGOoeTfBNX9RaSSKfk3Ydw4tfxFpHJpklcTtm4N\nY/43boTu3WNHIyLSNE3ySlC3buHJXlriWUQqkZJ/MzTkU0QqlZJ/M1T3F5FKpZp/MzZsgBEjwmJv\nVbpNikhKqeafsH79QqfvsmWxIxERSZaSfwtU+hGRSqTk3wJ1+opIJVLyb8H48fDss7GjEBFJlpJ/\nC44+GrZsgbq62JGIiCRHyb8FZjBhAjzzTOxIRESSk0jyN7PJZrbczFaY2XWN/PwrZrYwvz1nZscm\ncd1SUfIXkUpTdPI3syrgdmASMBKYZmYjGhy2Cpjg7scB/wbcVex1S+m005T8RaSyJNHyHwOsdPda\nd98OzAamFB7g7i+4+wf53ReAgQlct2SOPRbeeQfWr48diYhIMpJI/gOBwu7QNTSf3C8FHkvguiVT\nVQWnnqrWv4hUjs6lvJiZnQ5cAoxv7riZM2fufp/L5cjlch0aV2vUl36mTo0diYhkXXV1NdXV1UWd\no+i1fcxsLDDT3Sfn968H3N1vanDcKGAOMNnd32zmfKlZ26fQq6/ChRdCTU3sSERE9hZrbZ/5wHAz\nG2JmXYCpwCMNAhtMSPxfbS7xp9lxx4WF3lT3F5FKUHTyd/edwHRgLrAUmO3uNWZ2uZldlj/sB0Bv\n4A4zW2BmZfeIlE6dIJeDp56KHYmISPG0pHMbzJoVyj/33BM7EhGRPbSkcwebOBH+9CdI6b1JRKTV\nlPzbYMQI+PRTeOut2JGIiBRHyb8NzPa0/kVEypmSfxudcYaSv4iUP3X4tlFtLZx4Irz9tp7rKyLp\noA7fEhgyBHr1goULY0ciItJ+Sv7tMGkSPPFE7ChERNpPyb8dlPxFpNyp5t8OH30Ehx4alnro0SN2\nNCKSdar5l0j37jBmDDz9dOxIRETaR8m/nVT6EZFypuTfTkr+IlLOlPzbadSoUPtfsSJ2JCIibafk\n305mcPbZ8Ic/xI5ERKTtlPyLcPbZ8Pvfx45CRKTtNNSzCB99BP37Q10d9OwZOxoRySoN9Syx7t3h\n1FPV8Ssi5UfJv0gq/YhIOVLZp0h1dTB6dFjls1On2NGISBap7BPBYYeF7bnnYkciItJ6Sv4JOO88\neOih2FGIiLSeyj4JWLYszPitrdUDXkSk9FT2ieRznwure778cuxIRERaR8k/IX//9zBnTuwoRERa\nR8k/IfXJv8wqViKSUUr+CTn+eNi1S8/2FZHyoOSfEDP4x3+EBx6IHYmISMs02idBS5bAWWfBX/+q\nUT+l8uGH8NZbYVu1Kky2e+892Lw5vG7ZAjt2wM6d4dUdunXbs+2/P/TqBX36hK1v3z3v+/cPczgO\nOij271Kkee0Z7dM5oQtPBm4lfJO4x91vauSY24AzgY+Ai939tSSunSbHHAMHHxwmfE2YEDuayrNp\nE8yfH0ZV1W/vvw9Dh8KwYeG1f//wvnfvsPXoAZ07h61Tp/ANbds22LoVPvkEPv443CjefTdsixeH\n13feCc9orqsLv27QoHAjqH+t34YMgcGDw41EpJwU3fI3sypgBXAGsA6YD0x19+UFx5wJTHf3vzOz\nk4Afu/vYJs5Xti1/gBtvDC3/n/0sdiTlb/t2mDcPHn8cHnsstO5PPBFOOCFsX/hCSL7WpvZO27iH\nG8yaNeFGUP9av9XWhtdevUIsQ4bA4YfveV+/6duDdKT2tPyTSP5jgRnufmZ+/3rAC1v/ZvYz4Gl3\n/3V+vwbIufuGRs5X1sm/tjYkpXXroEuX2NGUn1274Nln4Re/gN/+NrTmzzwTJk+GsWNDCz5tdu0K\n3xJqa/fe/vrXPe+7dGn+5tC3b8fexKSyxSr7DATqCvbXAGNaOGZt/rN9kn+5GzIEjj4aHn0Uvvzl\n2NGUj9pauOsu+OUvQyv5oovgX/8VBg6MHVnLqqpCnAMHwimn7Ptz91BKanhT+POf93y2bVsoHw0a\nBP36Nb4dckjY9tuv5L/FzHHfe9h2Y+3Rhp+1tA/pahCmsB1V/i65BO69V8m/NV55BW65JTwT4aKL\n4OGH4bjjYkeVLLPQsu/bN3wrbMzf/garV4ey0saNsGFD2BYv3nv/nXega9fQt1S49ey553337qEj\nu6mta9fQj1FVFV6be79rV9h27mz9644d6dq2b2/7r3Hf95tYY9/MWjqmcP+II6Cmpu3/fjpKEsl/\nLTC4YH9Q/rOGxxzWwjG7zZw5c/f7XC5HLpcrNsaSuuAC+Pa3Q+lnwIDY0aTTn/4E//7v8MYbcPXV\noY8ky3Xxgw4KAwaOOab549zDCKb334cPPgivDbcPPww3jE8+aXzbtm1Poi5M2g3f79y550bQmtf6\n9/vtt6eTvZjtgAOSOU/h1trY0j5ar7q6murq6qLOkUTNvxPwOqHDdz3wEjDN3WsKjjkLuCLf4TsW\nuLVSO3zr/fM/h1En3/1u7EjS5ZVX4PrrQ6ljxoxwo1QZQ6Q4URZ2c/edwHRgLrAUmO3uNWZ2uZld\nlj/mUeAtM3sD+H/AN4u9btpdeincc09oSQm8+SZMnQrnnBOWwli6FC68UIlfJBZN8uog7jBqFNx2\nG5x+euxo4tm2DW66Kfw5XH01XHNNqEmLSHK0pHOKmME3vgG33x47knieeircABcsgFdfhe9/X4lf\nJC3U8u9AW7aEoZ+vvBLGdmfF5s1w1VXwzDPwk5/AuefGjkiksqnlnzI9esDFF8OsWbEjKZ25c0Nr\nv2fPUNdX4hdJJ7X8O1j9kgS1tZVd8vjoI7juujBO/9574Utfih2RSHao5Z9CQ4fCqafCfffFjqTj\nLFwIo0eHceeLFinxi5QDtfxLYN48mDYNVq6svKGN990H//IvcOutYeimiJRelIXdklaJyR9Ca3ja\nNPj612NHkoytW/d06s6ZAyNHxo5IJLtU9kmxH/wAfvjDsG5IuXvrLRg3LjwsZf58JX6RcqTkXyIT\nJoRVH8v9MY+PPhqWVv7qV+HBB+HAA2NHJCLtobJPCT3zTBj6WVNTfk9+2rkTbrghjOSZPRvGj48d\nkYjUU9kn5U47LSxXfNttsSNpm02bwgNVnn02TFhT4hcpf0r+JXbzzWHbtCl2JK3zwgthDfrRo+HJ\nJ8NDRUSk/KnsE8G3vhUeMPHTn8aOpGnucMcdodRz110wZUrsiESkKRrqWSY2bw4P7fj1r9NZQtmy\nBS67DJYtC8M4jzgidkQi0hzV/MtEr15hwbNLLw3j5dNk+XI46aTwqL9585T4RSqVkn8k550XWv83\n3BA7kj0efDAsRXHNNWFUz/77x45IRDqKyj4RbdgQOlLvvRcmTYoXxyefwLXXhhU5H3yw6YeMi0g6\nqexTZvr1g/vvh699DVavjhPD8uVh0tbmzeGBK0r8Itmg5B/ZaaeFVvc//EPoaC0Vd/jP/wxlnunT\nw8zjnj1Ld30RiUtlnxRwD52/dXXwP/8TOls70po14RGTq1fDr34Fxx7bsdcTkY6lsk+ZMoM774SD\nDgorf27b1jHXcQ9j9j//eRgzBl5+WYlfJKvU8k+RbdvCmvgbN8Jvfwt9+iR37mXL4Mor4cMPQwfz\nMcckd24RiUst/zLXtWsYbXPyyWGs/bx5xZ9z3Tq44grI5eCcc+D555X4RUTJP3WqquCmm+BHPwpz\nAa66Ct59t+3nqamBb34zJPpu3cL+1VdD587Jxywi5UfJP6XOPx+WLAkzgIcPDyNy5s0LSys35Y03\nwszh8eNh4kTo3TsM5bzllmRLSCJS/lTzLwPr14dF1h55JIzUGTkSBg8OLfqtW2Ht2lDT79QJJk8O\ni7CddVblPS9YRBqnhd0yYP16eP11qK2FTz8NN4BDDw03hP79w8ghEckWJX8RkQwq+WgfM+tlZnPN\n7HUze8LM9pkjamaDzOwpM1tqZovN7FvFXFNERIpXbIfv9cAf3f0o4Cngu40cswO41t1HAicDV5jZ\niCKvG1V1dXXsEFpFcSZLcSZLccZVbPKfAtyXf38f8OWGB7j72+7+Wv79FqAGGFjkdaMql38MijNZ\nijNZijOuYpP/Ie6+AUKSBw5p7mAzOxw4HnixyOuKiEgRWpzyY2ZPAoWP7TbAge83cniTPbVm1gP4\nDXBV/huAiIhEUtRoHzOrAXLuvsHMDgWedvejGzmuM/B74DF3/3EL59RQHxGRNmrraJ9iJ/s/AlwM\n3AR8DXi4iePuBZa1lPih7b8BERFpu2Jb/r2BB4HDgFrgAnd/38z6A3e5+9lmNg74M7CYUBZy4Hvu\n/njR0YuISLukbpKXiIh0vNQs7GZmk81suZmtMLPrYsfTmHKbsGZmVWb2qpk9EjuWpphZTzP7bzOr\nyf+5nhQ7pobM7BozW2Jmi8zsV2bWJXZM9czsHjPbYGaLCj5rcfJlCmK8Of93/pqZzTGzg2LGmI9p\nnzgLfvZtM9uVr3ZE1VScZnZl/s90sZnd2NJ5UpH8zawKuB2YBIwEpqV0Ili5TVi7ClgWO4gW/Bh4\nND9Q4DjCPJDUMLMBwJXAaHcfRegnmxo3qr38nPD/plBrJl+WUmMxzgVGuvvxwErixwiNx4mZDQK+\nRChtp8E+cZpZDjgHONbdjwX+o6WTpCL5A2OAle5e6+7bgdmECWSpUk4T1vL/YM8C7o4dS1Pyrb1T\n3f3nAO6+w93/FjmsxnQCuudHrR0ArIscz27u/hywucHHLU6+LKXGYnT3P7r7rvzuC8CgkgfWQBN/\nlgD/F/jfJQ6nSU3E+b+AG919R/6Yd1o6T1qS/0CgrmB/DSlNqvXKYMJa/T/YNHfqDAXeMbOf58tT\nd5rZ/rGDKuTu64BbgNXAWuB9d/9j3Kha1KbJlynwdeCx2EE0xszOBercfXHsWFpwJDDBzF4ws6fN\n7ISWfkFakn9ZSfuENTP7O2BD/luK5bc06gyMBma5+2jgY0LJIjXM7GBCS3oIMADoYWZfiRtVm6W2\nAWBm/wfY7u73x46loXxD5HvAjMKPI4XTks5AL3cfC3yHMAqzWWlJ/muBwQX7g/KfpU7+q/9vgP9y\n96bmNcQ2DjjXzFYBDwCnm9kvIsfUmDWEVtXL+f3fEG4GafJFYJW7v+fuO4GHgFMix9SSDWbWDyA/\n+XJj5HgaZWYXE0qTab2ZHgEcDiw0s7cIeekVM0vjN6k6wr9N3H0+sMvMmn1+X1qS/3xguJkNyY+k\nmEqYQJZGrZ6wFou7f8/dB7v7MMKf5VPuflHsuBrKlybqzOzI/EdnkL4O6tXAWDPrZmZGiDFVndLs\n++2ufvIlND/5spT2itHMJhPKkue6+7ZoUe1rd5zuvsTdD3X3Ye4+lNBY+by7p+Fm2vDv/HfARID8\n/6f93L3Zp3+nIvnnW1TTCSMAlgKz3T1t/8HIT1i7EJhoZgvyderJseMqc98CfmVmrxFG+/wwcjx7\ncfeXCN9IFgALCf/h7owaVAEzux94HjjSzFab2SXAjcCXzOx1ws2qxWF/EWL8CdADeDL//+iOmDFC\nk3EWclJQ9mkiznuBYWa2GLgfaLGxp0leIiIZlIqWv4iIlJaSv4hIBin5i4hkkJK/iEgGKfmLiGSQ\nkr+ISAYp+YuIZJCSv4hIBv1/sts9OwPgMpUAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.plot(ts, position, 'b')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Turning dimensions back on" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def wrapped_get_dimensionality_method(self):\n", " print('hi')\n", " return saved_get_dimensionality_method(self)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can restore the saved version of `UNITS._get_dimensionality`." ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ ">" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "UNITS._get_dimensionality = saved_get_dimensionality_method\n", "UNITS._get_dimensionality" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now the units and the quantities have dimensions again." ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(,\n", " ,\n", " ,\n", " )" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "UNITS.kilogram.dimensionality, UNITS.second.dimensionality, UNITS.newton.dimensionality, UNITS.meter.dimensionality" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(,\n", " ,\n", " ,\n", " )" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mass.dimensionality, k.dimensionality, c.dimensionality, ts.dimensionality" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`slope_func` still works." ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(, )" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "slope_func(y_init, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But `odeint` doesn't." ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false }, "outputs": [ { "ename": "ValueError", "evalue": "setting an array element with a sequence.", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0masol\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0modeint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mslope_func\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0my_init\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[1;32m/home/downey/anaconda2/lib/python2.7/site-packages/scipy/integrate/odepack.pyc\u001b[0m in \u001b[0;36modeint\u001b[1;34m(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg)\u001b[0m\n\u001b[0;32m 213\u001b[0m output = _odepack.odeint(func, y0, t, args, Dfun, col_deriv, ml, mu,\n\u001b[0;32m 214\u001b[0m \u001b[0mfull_output\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrtol\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0matol\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtcrit\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mh0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhmax\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhmin\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 215\u001b[1;33m ixpr, mxstep, mxhnil, mxordn, mxords)\n\u001b[0m\u001b[0;32m 216\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0moutput\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m<\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 217\u001b[0m \u001b[0mwarning_msg\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_msgs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0moutput\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;34m\" Run with full_output = 1 to get quantitative information.\"\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;31mValueError\u001b[0m: setting an array element with a sequence." ] } ], "source": [ "asol = odeint(slope_func, y_init, ts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.11" } }, "nbformat": 4, "nbformat_minor": 0 }