{ "metadata": { "name": "", "signature": "sha256:25c610c11390c8e8b45ab01dd178e174c77ba77b2a05d3c90751680b130234d0" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ " Numerical Differentiation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some material comes from the [Computational Physics](http://www-personal.umich.edu/~mejn/cp/) book by Mark Newman at University of Michigan, materials developed by [Matt Moelter](http://physics.calpoly.edu/node/88) and [Jodi Christiansen](http://physics.calpoly.edu/node/70) for PHYS 202 at [Cal Poly](http://physics.calpoly.edu/), as well as some [NumPy tutorials](http://wiki.scipy.org/Tentative_NumPy_Tutorial)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Instructions:** Create a new directory called `Differentiation` with a notebook called `DifferentiationTour`. Give it a heading 1 cell title **Numerical Differentiation**. Read this page, typing in the code in the code cells and executing them as you go. \n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Do not copy/paste. \n", "\n", "Type the commands yourself to get the practice doing it. This will also slow you down so you can think about the commands and what they are doing as you type them.\n", "\n", "Save your notebook when you are done, then try the accompanying exercises.\n", "\n", "---" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%pylab inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Introduction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Derivatives are everywhere in physics. Velocity is the time derivative of position. Acceleration is the time derivative of velocity. Derivatives tell us about the slope of a function and the minimum or maximum of a function occurs where the derivative is zero. We use these facts all the time so let\u2019s investigate some numerical derivatives and their applications." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You hear a lot less about numerical derivatives than integrals, for a number of reasons:\n", "\n", "1. The basic techniques for numerical derivatives are quite simple, so they don\u2019t take long to explain.\n", "\n", "2. Derivatives of known functions can always be calculated analytically, so there\u2019s less need to calculate them numerically.\n", "\n", "3. There are some significant practical problems with numerical derivatives, which means they are used less often then numerical integrals." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The standard definition of a derivative, the one you see in the calculus books, is\n", "\n", "$$\n", "\\frac{df}{dx} = \\lim\\limits_{h\\rightarrow 0} \\frac{f(x+h)\u2212f(x)}{h}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is the *forward difference* method - it\u2019s simply the slope of the curve $f(x)$ measured over a small interval of width $h$ in the forward direction from $x$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The backward difference, which has the mirror image definition, is given by\n", "\n", "$$\n", "\\frac{df}{dx} \\simeq \\frac{f(x) - f(x-h)}{h}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The forward and backward differences typically give about the same answer and in many cases you can use either. Most often one uses the forward difference. There are a few special cases where one is preferred over the other, particularly when there is a discontinuity in the derivative of the function at the point $x$ or when the domain of the function is bounded and you want the value of the derivative on the boundary, in which case only one or other of the two difference formulas will work. The rest of the time, you just pick whichever one is most convenient. A third option is the centered difference, which we'll learn about in a minute." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider a function like, $y = x^5$. Numerically in a computer, this function is represented as an array of $y$ values evaluated at $x$ values. " ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.arange(0.,10.01,0.1)\n", "y = x**5" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The numerical derivative is calculated using the pre-calculus method of finding $\\Delta y$ and dividing by $\\Delta x$. This is only a good approximation to the derivative if $\\Delta$ is small, but it\u2019s the best we\u2019ve got. We are just going to calculate the slope, $\\Delta y/\\Delta x$ at every point." ] }, { "cell_type": "code", "collapsed": false, "input": [ "dy = np.diff(y) #y(n+1)-y(n)\n", "\n", "#Check out the help on the diff method if you want to learn more:\n", "#help(np.diff)\n", "\n", "print len(y), len(dy)\n", "dx = np.diff(x)\n", "dydx = (dy/dx)\n", "print len(dydx)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following won\u2019t work. Can you see why?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "plt.plot(x,dydx)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The vector `x` is longer by one element than the derivative vector `dydx`. If we want to have the `dydx` with the same length as our original `x,y`, we can do:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "#specify the size of dy ahead because diff returns an array of n-1 elements\n", "dydx1 = np.zeros(y.shape,float) #we know it will be this size\n", "\n", "dydx1[0:-1] = np.diff(y)/np.diff(x)\n", "dydx1[-1] = (y[-1] - y[-2])/(x[-1] - x[-2])" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which will use the forward finite difference for all elements but the last one, where we use a backward difference. \n", "\n", "Recall the slicing syntax:\n", "\n", " array[:] means take the whole array\n", " array[2:7] means take the third through 6th elements (Python numbering)\n", " array[-1] means return the last element in the array (negative numbers count back from the end)\n", " array[1:-1] means take everything but the first and last elements in the array\n", "\n", "\n", "To appreciate how array slicing like `array[2:7]` works, think of the left-hand-side of the colon as one index and the right-hand-side as another and implement any logic for each side separately. So for example, to implement the two-point forward difference using any arbitrary start and end points, `i` and `j` you would use:\n", "\n", "```\n", "dydx[i:j] = (y[i+1:j+1] - y[i:j])/(x[i+1:j+1] - x[i:j])\n", "```\n", "\n", "As an alternative to the two-point forward difference, we can calculate the slope at the center of each x-bin. In this case, we have to treat the first and last elements specially, but we can use array slices for everything in between:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "#calculate dydx by center differencing using array slices\n", "dydx2 = np.zeros(y.shape,float) #we know it will be this size\n", "\n", "dydx2[1:-1] = (y[2:] - y[:-2])/(x[2:] - x[:-2]) #center difference\n", "\n", "dydx2[0] = (y[1]-y[0])/(x[1]-x[0]) #forward difference\n", "\n", "dydx2[-1] = (y[-1] - y[-2])/(x[-1] - x[-2]) #backward difference" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let\u2019s compare our numerical derivatives to the exact value: \n", " \n", "$$\\frac{dy}{dx} = 5x^4$$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "dydxExact = 5*x**4" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "plt.plot(x,dydx1,label='forward difference')\n", "plt.plot(x,dydx2,label='center difference')\n", "plt.plot(x,dydxExact,label='exact')\n", "plt.legend(loc='upper left')\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "They are *slightly* different but it is hard to see. Let\u2019s plot the error introduced by our numerical derivatives. Start from the second value in the array to avoid dividing by zero." ] }, { "cell_type": "code", "collapsed": false, "input": [ "percentError1 = 100*abs( dydx1[1:] - dydxExact[1:] )/dydxExact[1:]\n", "percentError2 = 100*abs( dydx2[1:] - dydxExact[1:] )/dydxExact[1:]\n", "\n", "#Use a semilog y-axis to see the variation better\n", "plt.semilogy(x[1:], percentError1, label='forward difference')\n", "plt.semilogy(x[1:], percentError2, label='center difference')\n", "plt.ylabel(\"percent error with dydxExact\")\n", "plt.legend(loc=\"upper right\")\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since the center difference is the more accurate of the two, we'll generally use that one. The overall accuracy of the derivative depends on the spacing of points. You can explore this phenomenon further with the accompanying exercises." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " All content is under a modified MIT License, and can be freely used and adapted. See the full license text [here](../../LICENSE)." ] } ], "metadata": {} } ] }