{ "metadata": { "name": "", "signature": "sha256:4e00fe6dbcc7a87ee494c321cd344e77e58d7baaff53f19dcc244a015f8c2bb8" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ " Interpolation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This lesson is based on 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)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "**Instructions:** Create a new directory called `Interpolation` with a notebook called `InterpolationTour`. Give it a heading 1 cell title **Interpolation**. Read this page, typing in the code in the code cells and executing them as you go." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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": [ "In this tour you will explore the concept of interpolation and learn how to use SciPy's built-in interpolation routines to approximate functions. A common use for interpolation is to smooth out a dataset that has relatively few datapoints." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Interpolation is a method of curve fitting a set of data to obtain values at unknown locations between two known points. There are different methods for interpolation, including the simplest - linear - which assumes that the function that connects the two known points is a straight line. Any values of the function for points between can be calculated from the interpolating function. Linear interpolation doesn't always give the best results, so higher order polynomial fits can be employed. Cubic is a common choice. \n", "\n", "For more on this topic, read up on it via [wikipedia](http://en.wikipedia.org/wiki/Interpolation).\n", "\n", "Before we use interpolation, let's review using boolean masks. We'll set up an example for finding data with a mask and then see how interpolation can help us when the data are sparse." ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Finding values with Boolean Masks" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sometimes we want to find a subset of the data points we\u2019ve computed. Here\u2019s an example. Lets look at the function $\\sin(\\theta)$ for $0\\lt \\theta\\lt2\\pi$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "theta = np.arange(0.,2*np.pi,np.pi/100)\n", "plt.plot(theta, np.sin(theta))\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now suppose that we want to plot just the data with $0\\lt \\sin{(\\theta)}\\lt 0.5$. We\u2019re going to have to \"find\u201d the data with $0\\lt \\sin{(\\theta)}$ and $\\sin{(\\theta)}\\lt 0.5$. We can use a mask to accomplish this. Set the mask to be the subset of the array that satisfies the criteria:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "ii = (np.sin(theta)>0) & (np.sin(theta)<0.5)\n", "plt.plot(theta,np.sin(theta),theta[ii],np.sin(theta[ii]),\"r-\")\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Oops, this has a horizontal line between the two sections of the `sin` function that pass our criteria. `ii` is now just a boolean mask where `theta` meets the logical criteria. Print the mask to see what it contains. It should be `True` in the array element locations passing our criteria and `False` everywhere else.\n", "\n", "Referring to `theta[ii]` will mean that we get just those elements of `theta` where the logical criteria was satisfied (where `ii` is `True`).\n", "\n", "Now we can get rid of the ugly line. `plot()` has decided to connect the 2 parts of the selected values with a straight line. To get two separate line segments we will need to set one value along the ugly line to `NaN` (\u201cNot a Number\u201d). For this case we will find the offending element where \u201c`y equals max(y)`\u201d and then set that element to `NaN`." ] }, { "cell_type": "code", "collapsed": false, "input": [ "y = np.sin(theta[ii])\n", "jj = (y == max(y))\n", "y[jj] = NaN\n", "plt.plot(theta, np.sin(theta))\n", "plt.plot(theta[ii],y,linewidth=4)\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Interpolation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What if our range is so small, that our boolean mask doesn\u2019t find anything? \n", "\n", "Create a new `theta` array that doesn't have a value of 0.5. When you try to find it with the mask, all entries in the mask will be `False`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "theta = np.arange(0.,2*np.pi+0.1,np.pi/10)\n", "ii = (theta == 0.5)\n", "print ii" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How can we numerically find a new value that is between the existing values? We\u2019ll use *interpolation*. \n", "\n", "We have to import the interpolation routines from `scipy` - the scientific python computing package." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.interpolate import interp1d \n", "#that's a \"one\" as in \"one-dimensional\", not an \"ell\"\n", "\n", "#Create the \"interpolating function\":\n", "sinInterp = interp1d(theta, sin(theta))\n", "\n", "#Set the x value where we want to know the function and \n", "#evaluate the interpolating function at that point:\n", "newTheta = 0.5\n", "newSin = sinInterp(newTheta)\n", "\n", "#Alternatively, here's a way to do it in one line:\n", "#newSin = interp1d(theta, np.sin(theta))(newTheta)\n", "\n", "#Now plot the interpolated point as a big plus sign:\n", "plt.plot(theta, np.sin(theta),newTheta, newSin, 'r+',markersize=20,markeredgewidth=3)\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "#Learn more about interp1d:\n", "help(interp1d)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could also find a whole bunch of new values if we use an array:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "thetaNew = np.arange(0,2*np.pi+0.1,np.pi/30)\n", "sinNew = interp1d(theta, np.sin(theta))(thetaNew)\n", "plt.plot(theta, np.sin(theta), thetaNew, sinNew, 'r.')\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let\u2019s look at how accurate our interpolated values are. For this we will compare the interpolated value to the exact value:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "#sinNew is the interpolated value, np.sin(thetaNew) is the exact value\n", "plt.plot(thetaNew, sinNew - np.sin(thetaNew))\n", "plt.ylabel('difference between interp and exact value')\n", "plt.xlabel(r'$\\theta$ (rad)',fontsize=20)\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Good to about 1/100. But it\u2019s easy to do better by using the curvature of the function instead of just straight line estimates between points. We\u2019ll ask for `cubic` interpolation. Review the [documentation](http://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html) to learn more about how this works." ] }, { "cell_type": "code", "collapsed": false, "input": [ "#cubic interpolation\n", "sinNewC = interp1d(theta, np.sin(theta), kind='cubic')(thetaNew)\n", "\n", "#Plot them both together\n", "plt.subplot(2,1,1)\n", "plt.plot(thetaNew, sinNew - np.sin(thetaNew),thetaNew, sinNewC - np.sin(thetaNew), 'r')\n", "\n", "#Can't really see the variation in the red curve. Let's zoom in:\n", "plt.subplot(2,1,2)\n", "plt.plot(thetaNew, sinNewC - np.sin(thetaNew),'r')\n", "\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Yowza, that is quite a lot better." ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Another example" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.arange(0, 10) \n", "y = np.array([3.0, -4.0, -2.0, -1.0, 3.0, 6.0, 10.0, 8.0, 12.0, 20.0]) \n", "f = interp1d(x, y, kind = 'cubic') \n", " \n", "xint = 3.5 #location to evaluate\n", "yint = f(xint) #interpolated point at that location\n", " \n", "plt.plot(x, y, 'o', c = 'b') \n", "plt.plot(xint, yint, 's', c = 'r') \n", " \n", "plt.show() \n" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.arange(0, 10) \n", "y = np.array([3.0, -4.0, -2.0, -1.0, 3.0, 6.0, 10.0, 8.0, 12.0, 20.0]) \n", "flin = interp1d(x, y) #linear is the default\n", "fcub = interp1d(x, y, kind = 'cubic')\n", "\n", "#Create a full interpolated curve\n", "xint = np.arange(0, 9.01, 0.01) \n", "ylin = flin(xint) \n", "ycub = fcub(xint)\n", "\n", "plt.subplot(211)\n", "plt.plot(x, y, 'o', c = 'b') \n", "plt.plot(xint, ycub, '-r',label=\"cubic\") \n", "plt.legend(loc=\"best\")\n", "plt.subplot(212)\n", "plt.plot(x, y, 'o', c = 'b') \n", "plt.plot(xint, ylin, '-g',label=\"linear\") \n", "plt.legend(loc=\"best\")\n", "\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Higher dimensions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Boolean masks work in 2D too. Recall the point potential function from the Graphics tour. Let's put it into a function that we can re-use. We'll save it to a file called `Electrostatics.py` and import it as a module:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%%file Electrostatics.py\n", "\n", "import numpy as np\n", "\n", "def point_potential(x,y,q,Xc,Yc):\n", " \"\"\"\n", " Return the electric potential at (x,y)\n", " for a point charge q at (Xc,Yc)\n", " \n", " Units returned are [Volts] if input \n", " units are [meters] and [Coulombs] \n", " \"\"\"\n", " k = 8.987551787997912e9 #(Nm^2/C^2) \n", " Vxy = k*q/np.sqrt(((x-Xc)**2 + (y-Yc)**2))\n", " return Vxy\n", "\n", "def dipole_potential(x,y,q,d):\n", " \"\"\"\n", " Return the electric potential at (x,y) for a \n", " pair of point charges +/- q \n", " separated by distance d along the x-axis \n", " with their midpoint at (0,0). \n", " \n", " Units returned are [Volts] if input\n", " units are [meters] and [Coulombs] \n", " \"\"\"\n", " Vxy=point_potential(x,y,q,-d/2.,0.) + point_potential(x,y,-q,+d/2.,0.)\n", " return Vxy\n", "\n" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "from Electrostatics import *" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "#Check out the help we created with those docstrings:\n", "help(point_potential)\n", "help(dipole_potential)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "#Create a couple of point potentials\n", "x = np.arange(-5.,5.01,0.51)\n", "y = np.arange(-5.,5.01,0.51)\n", "xg,yg = np.meshgrid(x,y)\n", "V1 = point_potential(xg,yg,1e-9,2.01,-2.01);\n", "V2 = point_potential(xg,yg,-1e-9,-2.01,2.01);" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "#Now plot in 3D\n", "from mpl_toolkits.mplot3d import Axes3D\n", "fig = plt.figure()\n", "ax = fig.gca(projection='3d')\n", "ax.plot_surface(xg,yg,V1+V2,rstride = 1, cstride = 1,cmap=cm.coolwarm);" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To find where `V1+V2 = 0`, we can find where `(V1 == -V2)`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "#Make the mask\n", "ii = (V1 == -V2)\n", "\n", "#Plot the x,y values using the mask:\n", "plt.plot(xg[ii], yg[ii]) \n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These are the $(x,y)$ values where the potentials exactly cancel." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Of course, we could also make a contour plot and see the contour at `V1+V2 = 0`, but this method with the mask gives us access to the array elements for use in further calculations." ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "`Interp2d`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2D interpolation is also similar to the 1D analog. A simple way to increase the sampling is to make a new mesh of points and then use `interp2d` to get the 2D interpolation:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.interpolate import interp2d\n", "xnew = np.arange(-5.,5.01,0.1)\n", "ynew = np.arange(-5.,5.01,0.1)\n", "xng,yng = np.meshgrid(xnew,ynew)\n", "#V1+V2 is the old grid using x,y\n", "Vnew = interp2d(x,y,V1+V2,kind='linear')(xnew, ynew)\n", "\n", "fig = figure()\n", "ax = fig.gca(projection='3d')\n", "ax.plot_surface(xng,yng,Vnew,cmap=cm.coolwarm);" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You\u2019ll notice that this procedure doesn\u2019t reproduce peaks very well, but in between the peaks it's not too bad. Let\u2019s check the accuracy of the calculation." ] }, { "cell_type": "code", "collapsed": false, "input": [ "fig = plt.figure()\n", "ax = fig.gca(projection='3d')\n", "ax.plot_surface(xng,yng,Vnew-point_potential(xng,yng,1e-9,2.01,-2.01) \\\n", " -point_potential(xng,yng,-1e-9,-2.01,2.01));" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is pretty flat except right near the poles." ] }, { "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": {} } ] }