gistfile1.txt

{
 "metadata": {
  "name": "",
  "signature": "sha256:3c570ceaf2bc0daef62dbd3b9dc80cf44ed37b93859610e9f5c30fec19d563ce"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Plot data points and functions"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline\n",
      "x = [1, 2, 3]\n",
      "y = [6, 11, 18]\n",
      "plt.plot(x,y,'bo',markersize=10)\n",
      "plt.xlim(0, 4)\n",
      "plt.ylim(5, 19)\n",
      "x=np.linspace(0, 4, 20)\n",
      "y=x**2+2*x+3\n",
      "plt.plot(x, y, 'r--')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "[<matplotlib.lines.Line2D at 0x106892b10>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAD7CAYAAACYLnSTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFmdJREFUeJzt3X10XHWdx/F3bdPSLKCCbMGWblYtuotYQVMahHSssAcs\nBxBjXRaky5MPoCIgCrJC6gIKFTiLCisPLfYI1RIQ2vL8NARsWgJIeBSsh5IallbpA9Spbdpk/7jT\nNoakM5mn352Z9+ucnM7cuXPnc35tv7nzu997L0iSJEmSJEmSJEmSJEmSJBXMsGJteOLEib0dHR3F\n2rwkVaoO4ONDfdO7ihAEgI6ODnp7e2P/c9FFFwXPYE5zlmtGcxb+B5iYS80tWjGXJJWOxVxS8cyf\nD+vXh05RFaq+mCcSidARsmLOwiqHnOWQEXaQ87HH4KyzYPPmkuYZTLmMZ66KdgAU6E3P/0iqNhs2\nwMSJcPnlcMwxodOUlWHDhkEOtdliLqnwzj8fli2DW28NnaTs5FrMRxQ+iqSq9vTTcOON8OyzoZNU\nlaqfM5dUYHfdBbNmwZ57hk5SVZxmkaQYyXWaxT1zSaoAmebMZwPTgFXAfullk4CfAjXAZuB0oL1Y\nASXFTyqVorV1CS0tbXR2bmD8+NE0NTXQ2DiZ2tra0PGqUqZd+UOA9cBcthfzJPBD4D7gCOA7wKcH\neK/TLFIFmj79QtrbR9HV1UB392SgFkhRU7OEsWPbqK/fyPz5Pwgds2wVq5vlMaCu37L/A96dfvwe\noGuoHyqpPKVSKdrbR7F8+QXblo1jBbuxmme7p7J8+VTgYlKplHvoJZbLnPl5wBVAJzALOL+giSTF\nVmvrErq6Gvos6eXnfIWjWLBtSVdXA62tS0ofrsrlUsxvBL4JjAfOIppXl1QFWlra0lMrkf/gFsbS\nxWV8d9uy7u4GWlraQsSrarmcNDQJODT9uAW4YbAVm5ubtz1OJBIVf20EqdJ1dm4gmiOHPVjFlZzN\nNO6im5F91qpNr6dsJJNJkslk3tvJpZgvA6YAjwJTgVcGW7FvMZdU/saPHw2kgFr+hzOZy4k8xSf7\nrZVKr6ds9N/RnTlzZk7byVTM5xEV7vcBK4ALgS8DPwNGARvSzyVVgaamBubOXcIe3R9mAn/gFG58\nxzo1NW00NTUM8G4Vk2eASspaKpVi332vYvnyCxhGD70DHHarq7uYF144226WHHmhLUlFV1tbS339\nRuDidJ95A9v7zNvSfeabLOQBuGcuacg8A7R4vJ65JFUAL7QlqbiSSXAHLbYs5pIyu/demDHDmzPH\nmAdAJe3YmjVw6qkwdy7sskvoNBqEc+aSduz442H33eHqq0MnqQq2JkoqvJYWePJJ+N3vQidRBs6Z\nSxrc3Lnwi1+A7Yax5zSLpMH19sKwYpYJ9WdroqTCs5CXDYu5JFUAi7kkVQCLuaTt7r0X/vKX0CmU\nA4u5pMjLL8MJJ0QnCansWMwlwebNcOKJ0NwMEyaETqMcWMwlweWXR6fqn3566CTKkX3mUrXr6IBD\nD4Wnn4a99w6dpurZZy4pN488AldcYSEvc+6ZS1KMuGcuSVXMYi5JFcBiLkkVIFMxnw2sBJ7rt/wb\nwEvA88BlRcglqVhaW2Hx4tApVGCZivkc4PB+yz4NHAV8DPgo8OMi5JJUDGvXRmd5vv126CQqsGyO\nmNYBC4H90s/nA/8LPJzhfXazSHEzY0Z0o4lrrw2dRIMo5W3jJgCNwKXA34BvA0/msB1JpXTHHfDb\n38Izz4ROoiLIpZiPAN4LTAbqifbUPzDQis3NzdseJxIJEolEDh8nKW+rVsHXvga33go77xw6jfpI\nJpMkk8m8t5PLNMs9wI+AR9PPlwEHAm/2e5/TLFJctLXBo4/CeeeFTqIMSjnNcgcwlaiY7wOM5J2F\nXFKcNDREP6pYmYr5PGAKsDuwAriQqF1xNlG74ibgxGIGlCRl5rVZJClGvDaLJFUxi7lUiW6/Ha67\nLnQKlZDFXKo0nZ3w1a/CxImhk6iELOZSJdmyJTpd/6yz4MADQ6dRCVnMpUpyySUwYgR85zuhk6jE\ncukzlxRHv/0tXHMNPPUUDB8eOo1KzNZEqVKsWAGvvgqNjaGTKA+5tiZazCUpRuwzl6QqZjGXpApg\nMZfKVXc39PSETqGYsJhL5er88+HSS0OnUEzYmiiVo/vug1//2rsGaRuLuVRuVq2Ck06Cm2+G3XcP\nnUYxYWuiVE56euDII2H//aOzPVVxbE2UqsG8ebB6NfS5v64E7plL5WXzZlizBvbYI3QSFYlngEpS\nBXCaRZKqmMVckiqAxVyKs9degz//OXQKlQGLuRRXGzbAMcfAb34TOonKgAdApbg67TR4++2oHXFY\nMf+rKk6KdQB0NrASeG6A184BeoDdhvqhkjK46SZ4/HG4/noLubKSqZjPAQ4fYPnewGHAawVPJFW7\nZ5+Fc8+FlhbYZZfQaVQmMhXzx4A1Ayy/EvCOsVIxPPIIXHUV7Ltv6CQqI7lcaOto4E/AswXOIgng\nzDNDJ1AZGmoxrwW+RzTFstWgE3rNfa4fkUgkSCQSQ/w4SapsyWSSZDKZ93ayObJSBywE9kv/PAik\n0q+NA7qAScCqfu+zm0WShijXbpah7pk/B4zp8/xV4BPA6qF+sCSpcDIdAJ0HLAb2AVYAJ/V73V1v\nKR9btsCXvgR//GPoJCpzmfbMj8vw+gcKFUSqSs3N8PrrUFcXOonKnLeNk0K55x6YMweeegqGDw+d\nRmXO0/mlEF57DSZNgttug4MPDp1GMeLNKaRy0dMDBx0EX/gCnHNO6DSKGYu5VE7a2+GTn/S6K3oH\ni7kkVQBvGydJVcxiLkkVwGIuFVtvL7zwQugUqnAWc6nYrr0WZsyIulikIvEAqFRMTzwBRx4JixfD\nhz4UOo3KgAdApbh5802YPh2uu85CrqJzz1wqhs2b4YgjYP/94fLLQ6dRGXHPXIqTZctgjz3g0ktD\nJ1GVcM9ckmLEPXNJqmIWc0mqABZzqRCcUlRgFnMpX7290a3fFi0KnURVzGIu5euyy+Dll+Eznwmd\nRFXM28ZJ+Vi0CH7yE1i6FEaPDp1GVczWRClXL74IiQQsWACTJ4dOowrhzSmkUurthcZGOPXU6CJa\nUoFYzKVSW7cO3v3u0ClUYYp50tBsYCXwXJ9ls4CXgA7gdsB/0ao+FnLFSDbFfA5weL9l9wP7AhOB\nV4DzC5xLkjQE2RTzx4A1/ZY9AGy90v5SYFwhQ0mShqYQfeYnA3cXYDtSfHV1waxZoVNIg8q3z/wC\nYBNwy0AvNjc3b3ucSCRIJBJ5fpwUwIYNcMwx8PnPh06iCpRMJkkmk3lvJ9sjpnXAQmC/Psv+EzgN\n+AzwtwHeYzeLyl9vL5xwQvTnzTfDsGI2gEm5d7Pkumd+OHAuMIWBC7lUGWbNik7Vf+wxC7liLZt/\nnfOIivb7iFoULyLqXhkJrE6v0wac3u997pmrvN13H5x8cnSq/jiP8as0PGlIKrSVK6MDnwccEDqJ\nqojFXJIqgLeNk6QqZjGXpApgMZe2WrAAtmwJnULKicVcArjmGjj3XHjrrdBJpJx4AFRauBC+/GV4\n/HH44AdDp1GVK/VJQ1LZSaVStLYuoaWljc7ODYwfP5qTP7YbDf99IcMWLbKQq6y5Z66qMH36hbS3\nj6Krq4Hu7slALf/ESyzmEJr/8XDWTvkA8+f/IHRMydZEaTCpVIr29lEsX34B3d1TgVoA1rEn3+Dn\nXL/ql7S3jySVSoUNKuXBYq6K19q6hK6uhncsX8t7uZ3oSohdXQ20ti4pdTSpYCzmqngtLW3pqZXB\ndXc30NLSVqJEUuFZzFXxOjs3sHVqZXC16fWk8mQxV8UbP340kGIqD1HDpkHWSqXXk8qTxVwVr6mp\ngWOGX8kvOYE9eWPAdWpq2mhqeue8ulQu7DNXxUvsPJL63h/yWR5mBeMHXGfs2DYaG88ucTKpcCzm\nqmyvvspO06dz9eQjWfX6A9R0rae7u4FoDj1FTU0bY8e2UV+/idraTPPqUnx50pAq1+rV8KlPwRln\nwNe/PuAZoE1NDTQ2TraQKza8OYXU3/r1MH9+dOs3qUxYzCWpAng6vyRVMYu5JFUAi7kqxyOPwAbP\n4lR1spirMixcCMcdBytWhE4iBZGpmM8GVgLP9Vm2G/AA8ApwP/Ce4kSTsvTww3DKKVFB32ef0Gmk\nIDIV8znA4f2WnUdUzPcBHko/l8JYuhS++MWoBbG+PnQaKZhs2l/qgIXAfunnvwemEO2x7wkkgY8M\n8D5bE1Vcr7wChxwCs2fDtGmh00gFUcp7gI4hKuSk/xyTwzak/I0fH+2RT5kSOokUXL7XZulN/wyo\nubl52+NEIkEikcjz46Q+dtrJQq6yl0wmSSaTeW8n12mWBPAGsBfwCE6zSFJBlPIM0AXAjPTjGcAd\nOWxDklRAmYr5PGAx8GFgBXAS8CPgMKLWxKnp51Jx/fWv8O1ve1KQNIhMc+bHDbL80EIHkQa1cSMc\neyy8//0walToNFIsedVExdvmzVEf+bBh8KtfwQjvp6LKVsrWRKk0enrg1FOj65IvWGAhl3bA/x2K\nr9mzYdkyuO8+p1ekDJxmUXx1d0cHPHfdNXQSqWS805AkVQDvNCRJVcxirvjwm5yUM4u54uHOO6Ne\nckk5sZtF4d12G5x+Otx9d+gkUtlyz1xhzZ8PX/961H74iU+ETiOVLYu5wrnlFvjWt+D+++HjHw+d\nRiprTrMonCeegAcegH33DZ1EKnv2mUtSjNhnLklVzGIuSRXAYq7SuOUWWLs2dAqpYlnMVXxXXAHf\n/z689VboJFLFsptFxXXZZXDDDZBMwt57h04jVSyLuYrnkktg7tyokI8dGzqNVNEs5iqOO++Em2+O\nCvlee4VOI1U8+8xVHFu2wLp1sNtuoZNIZcWbU0hSBfCkIUmqYvkU8/OBF4DngFsA77hbrXp74e23\nQ6eQqlquxbwOOA04ANgPGA78e4EyqZz09MAZZ8BZZ4VOIlW1XLtZ3gK6gVpgS/rPrkKFUpnYuBFm\nzIDXX4dFi0Knkaparnvmq4ErgE7gdWAt8GChQqkMvPUWTJsGmzZFN5bYddfQiaSqluue+QeBbxFN\nt6wDbgWOB27uu1Jzc/O2x4lEgkQikePHKVbWrYNPfxomTYKf/QyGDw+dSCpbyWSSZDKZ93ZybU38\nInAYcGr6+ZeAycAZfdaxNbFS9fbCggVw1FEwrJjdrVL1KXVr4u+Jivfo9IceCryY47ZUboYNg6OP\ntpBLMZJrMe8A5gJPAs+ml11XkESSpCHzDFBltnEjjPI0AqkUPANUxXHlldHcuKRY86qJGlhPD5x3\nHixcGLUeSoo1i7neqbsbTjkFli2Dxx+H3XcPnUhSBhZz/b1Nm6JOlREj4MEHobY2dCJJWfAAqP5e\nby/8+tfQ1BQVdEkl5fXMJakC2M0iSVXMYl7tNm8OnUBSAVjMq9n8+XDwwVEboqSy5hGuatTTAzNn\nwk03RRfMepe/06VyZzGvNuvXw4knwsqV8MQTMGZM6ESSCsBdsmqSSsGnPgXvfS88/LCFXKogtiZW\nm8cfjwq6l6+VYsk+c0mqAPaZS1IVs5hXqjffhOefD51CUolYzCvRCy/AgQdGl6+VVBVsTaw0ixbB\nSSfBj38MM2aETiOpRCzmlaK3Fy6/HK6+Otojnzw5dCJJJWQxrxRPPAEtLbB0KYwbFzqNpBKzNbGS\nbNkCw4eHTiEpD7YmykIuVbF8plneA9wA7Av0AicDSwoRSpFUKkVr6xJaWtro7NzA+PGjaWpqoLFx\nMrXezk1SH/lMs/wCeBSYTfRL4R+AdX1ed5olD9OnX0h7+yi6uhro7p4M1AIpdh7RytWjf8Cy+o9w\nyUOzQ8eUVGClnmZ5N3AIUSEH2MzfF3LlIZVK0d4+iuXLL6C7eypRIYcP0cWjm7/Hrm/vxaI/7EUq\nlQobVFJs5FrM/xn4MzAHeBq4nq0VR3lrbV1CV1fD3y07jltYzEHcwKk00cJLb0yltdVZLUmRXIv5\nCOAA4Jr0n38FzitUqGrX0tKWnlqJzOLbNNPMv3E/13I6MIzu7gZaWtrChZQUK7keAP1T+qc9/byF\nAYp5c3PztseJRIJEIpHjx1WXzs4N9P2iM5/pzOQi1rNLn7Vq0+tJKmfJZJJkMpn3dnIt5m8AK4B9\ngFeAQ4EX+q/Ut5gre+PHjwZSbC3o7UwaYK1Uej1J5az/ju7MmTNz2k4+febfAG4GOoCPAZfmsS31\n0dTUQE3NjufDa2raaGpq2OE6kqpHPsW8A6gHJgLHYjdL/p55BubMobFxMmPH7ng+fOzYNhobvf6K\npIhngMZBby9cey0cdhjstBO1tbXU12+kru5iamoeIppyAUhRU/MQdXUXU1+/yROHJG3jtVlCW7sW\nTjsNli2D+fNhwoRtL3kGqFR9vAdoOerogM99DqZNg1mzYKedQieSFJjFvBy99ho89RQce2zoJJJi\nwmIuSRXAS+BKUhWzmJfCmjXRPTn9piKpSCzmxbZgAXz0o9DZCd3dodNIqlBVX8wLcU2EAb35Jhx/\nPJxzDsybF91oeeTInDdXtJwFZs7CKYeMYM64sJgX4y/4xRejvfExY6L2w8bGvDdZLv8QzVk45ZAR\nzBkX+dw2ToOZMCGaXqmvD51EUpWo+j3zoqipsZBLKqli9pk/Q3QRLklS9jqAj4cOIUmSJEmSpJwd\nDvwe+APw3UHWuTr9egewf4ly9ZcpZ4LoBhu/S//8V8mSbTcbWAk8t4N14jCWmXImCD+WewOPEN3O\n8Hngm4OsF3o8s8mZIPx47gQsJToW9iLww0HWCz2e2eRMEH48txqezrBwkNdLNp7DgWVAHVBDNID/\n0m+dzwJ3px8fCOz4fmjFkU3OBLCgpKne6RCiv7DBimQcxhIy50wQfiz3ZPtBpJ2Bl4nnv81sciYI\nP56w/S7jI4jG6uB+r8dhPCFzzgTxGE+As4luvzlQniGNZ76tiZOIiuRyoBv4FXB0v3WOAn6RfrwU\neA8wJs/PHapsckJxu3uy8RiwZgevx2EsIXNOCD+WbxD90gZYD7wEvL/fOnEYz2xyQvjxhO23vBpJ\ntIO0ut/rcRhPyJwT4jGe44gK9g0MnGdI45lvMR8LrOjz/E/pZZnWGZfn5w5VNjl7gYOIvs7cDfxr\naaINSRzGMhtxG8s6om8SS/stj9t41jFwzriM57uIfvGsJJoaerHf63EZz0w54zKeVwHnAj2DvD6k\n8cy3mGd7GcD+v3VKffnAbD7vaaL5y4nAT4A7ipood6HHMhtxGsudgRbgTKI93/7iMp47yhmX8ewh\nmhIaBzQSTVf0F4fxzJQzDuN5JLCKaL58R98Ssh7PfIt5F9GgbLU30W+PHa0zLr2slLLJ+Tbbv57d\nQzS3vlvxow1JHMYyG3EZyxrgNuCXDPwfNi7jmSlnXMZzq3XAXcAn+y2Py3huNVjOOIznQUTTKK8C\n84CpwNx+65R0PEcAfyT6ejiSzAdAJxPmoEg2Ocew/bfgJKL59RDqyO4AaKix3KqOwXPGYSyHEf3n\nuGoH68RhPLPJGYfxfB/RnC3AaKAV+Ey/deIwntnkjMN49jWFgbtZSj6eRxAdgV8GnJ9e9pX0z1Y/\nTb/eARxQ7ECDyJTzDKLWsGeAxUSDV2rzgNeBTURzZScTz7HMlDMOY3kw0dftZ9jegnYE8RvPbHLG\nYTz3I5qeeAZ4lmiuF+I3ntnkjMN49jWF7d0scRtPSZIkSZIkSZIkSZIkSZIkSZIkSaH8P3xiYlTY\nd2vGAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x105910750>"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Polynomial interpolation using Newton's method "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "\n",
      "def evalPoly(a,xData,x):\n",
      "    n = len(xData)-1  # order of polynomial\n",
      "    p = a[n]\n",
      "    for i in range(1,n+1):\n",
      "        #print \"a[\", (n-i),\"]=\", a[n-i]\n",
      "        p = a[n-i] + (x-xData[n-i])*p\n",
      "    return p\n",
      "\n",
      "def coefficients(xData,yData):\n",
      "    n = len(xData)\n",
      "    a = np.zeros((n, n));\n",
      "    for i in range(0, n):\n",
      "        a[i, 0] = yData[i];\n",
      "    for j in range(1, n):\n",
      "        for k in range(0, n-j):\n",
      "            a[k, j] = (a[k+1, j-1] - a[k, j-1])/(xData[k+j]-xData[k])\n",
      "    return a[0, :] # return the zeroth row\n",
      "\n",
      "# y = x**2 + 2*x + 3\n",
      "xData = [1, 2, 3]\n",
      "yData = [6, 11, 18]\n",
      "a = coefficients(xData, yData)\n",
      "print a\n",
      "print \" x    yInterp   yExact\"\n",
      "print \"-----------------------\"\n",
      "for x in np.arange(0.0,4,0.2):\n",
      "    y = evalPoly(a,xData,x)\n",
      "    yExact = y = x**2 + 2*x + 3\n",
      "    print \"%3.1f %9.5f %9.5f\"% (x,y,yExact)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 6.  5.  1.]\n",
        " x    yInterp   yExact\n",
        "-----------------------\n",
        "0.0   3.00000   3.00000\n",
        "0.2   3.44000   3.44000\n",
        "0.4   3.96000   3.96000\n",
        "0.6   4.56000   4.56000\n",
        "0.8   5.24000   5.24000\n",
        "1.0   6.00000   6.00000\n",
        "1.2   6.84000   6.84000\n",
        "1.4   7.76000   7.76000\n",
        "1.6   8.76000   8.76000\n",
        "1.8   9.84000   9.84000\n",
        "2.0  11.00000  11.00000\n",
        "2.2  12.24000  12.24000\n",
        "2.4  13.56000  13.56000\n",
        "2.6  14.96000  14.96000\n",
        "2.8  16.44000  16.44000\n",
        "3.0  18.00000  18.00000\n",
        "3.2  19.64000  19.64000\n",
        "3.4  21.36000  21.36000\n",
        "3.6  23.16000  23.16000\n",
        "3.8  25.04000  25.04000\n"
       ]
      }
     ],
     "prompt_number": 55
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example from the textbook: first order polynomial interpolation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline\n",
      "\n",
      "def evalPoly(a,xData,x):\n",
      "    n = len(xData)-1  # order of polynomial\n",
      "    p = a[n]\n",
      "    for i in range(1,n+1):\n",
      "        #print \"a[\", (n-i),\"]=\", a[n-i]\n",
      "        p = a[n-i] + (x-xData[n-i])*p\n",
      "    return p\n",
      "\n",
      "def coefficients(xData,yData):\n",
      "    n = len(xData)\n",
      "    a = np.zeros((n, n));\n",
      "    for i in range(0, n):\n",
      "        a[i, 0] = yData[i];\n",
      "    for j in range(1, n):\n",
      "        for k in range(0, n-j):\n",
      "            a[k, j] = (a[k+1, j-1] - a[k, j-1])/(xData[k+j]-xData[k])\n",
      "    return a[0, :] # return the zeroth row\n",
      "\n",
      "xAll = [0, 10, 15, 20, 22.5, 30]\n",
      "yAll = [0, 227.04, 362.78, 517.35, 602.97, 901.67]\n",
      "plt.plot(xAll,yAll,'bo',markersize=10)\n",
      "\n",
      "xData = xAll[2:4]\n",
      "print \"xData=\",xData\n",
      "yData = yAll[2:4]\n",
      "print \"yData=\",yData\n",
      "a = coefficients(xData, yData)\n",
      "print \"a=\",a\n",
      "y = evalPoly(a, xData, 15)\n",
      "print \"f(15)=\", y\n",
      "y = evalPoly(a, xData, 16)\n",
      "print \"f(16)=\", y\n",
      "\n",
      "xs=np.linspace(0, 30, 50)\n",
      "ys=np.zeros(50)\n",
      "i=0\n",
      "for x in xs:\n",
      "    ys[i] = evalPoly(a,xData,x)\n",
      "    i=i+1\n",
      "plt.plot(xs,ys,'r--')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "xData= [15, 20]\n",
        "yData= [362.78, 517.35]\n",
        "a= [ 362.78    30.914]\n",
        "f(15)= 362.78\n",
        "f(16)= 393.694\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 54,
       "text": [
        "[<matplotlib.lines.Line2D at 0x10c336610>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGo1JREFUeJzt3Xt4FPW9x/F3wFWJKNSDh8sCxgPIrYKIQaIWAlUazpGL\nGlGr1qO0T+ulajlHFFFIBVrFo8g5PhQFbbEVWhqtRVtRRNZYm0hQEQyXJOgaskoUFVCXSwh7/vhN\nzCYkIZvdzezMfF7Ps8/Ozs5uvvMMzHfn9/3N7wciIiIiIiIiIiIiIiIiIiIiIiIiIuJiTwFVwOao\ndacCa4BS4BWgc9R7M4AyYBswLmr9cOs7yoCFSYxXREQS5HvAMOongPnAdGv5LuABa3kQsBHwARlA\nOZBmvbceGGEt/x3ISVrEIiKSMBnUTwDbgK7WcjfrNZhf/3dFbbcaGAl0B7ZGrb8KWJyMQEVE5Nja\nxfHZrphmIazn2mTQA6iM2q4S8DeyPmStFxERG8STAKJFrIeIiDjEcXF8tgrT9LML07zzqbU+BPSK\n2q4n5pd/yFqOXh9q7Iv79OkT2bFjRxyhiYh40g6gb0s3jucKYBVwvbV8PfB81PqrgOOBM4B+mOLv\nLmAfcB6mKHxd1Gfq2bFjB5FIxLWP2bNn2x6D9k37p/1z3wPoE8tJvKVXACuA0UAXYCcwC9PrZyUw\nFQgCU6xtt1jrtwCHgZupax66Gfgd0AHTC2h1LMGKiEjitDQBXN3E+ouaWP8r69HQ28BZLfybIiKS\nRIkqAksMsrOz7Q4hady8b6D9czq371+s0o69iS0iVnuWiIi0UFpaGsRwXtcVgIiIRykBiIh4lBKA\niIhHKQGIiHiUEoCIiEcpAYiIeJQSgIiIR8UzGJyIiNgsHA5TUFBEfn5hzJ/VjWAiIg41ZcosiotP\nIBTKorp6JHASxHBe1xWAiIgDhcNhiotPIBic2ervUA1ARMSBCgqKCIWy4voOJQAREQfKzy+0mn1a\nTwlARMSBKir2A+lxfYcSgIiIA/Xu3QEIx/UdSgAiIg6Um5uFz1cU13coAYiIONCoUSPx+2Pv+x9N\nCUBExIHS09PJzDxIRsZcfL61tKY5SDeCiYg4WPSdwE8+eS/EcF5XAhARcQlNCSkiIi2iBCAi4lFK\nACIiHqUEICLiUUoAIiIepQQgIuJRmg9ARMTpPvwQSktj/lgirgBmACXAZmA5cAJwKrAGKAVeATo3\n2L4M2AaMS8DfFxHxpspKuOkmOPdcKCmJ+ePxJoAM4CfAOcBZQHvgKuBuTAI4E1hrvQYYBFxpPecA\nixIQg4iI9yxZAkOGwCmnwPbtMG1azF8RbxPQPqAaMyh1jfX8MeZX/mhrm2VAAJMEJgErrM8EgXJg\nBBDfkHYiIl4zZoz51d+9e6u/It5f318ADwMVmBP/Hswv/65AlbVNlfUaoAdQGfX5SsAfZwwiIt7T\nt29cJ3+IPwH0Ae7ANAX1ADoC1zbYJmI9mqJBf0REGvPNN/DggxAKJeXr420COhf4J/C59fo5IAvY\nBXSznrsDn1rvh4BeUZ/vaa07Sl5e3rfL2dnZZGdnxxmqiIhDHDgAixebk//o0XD11Y1uFggECAQC\nrf4z8Y4GOhR4BsgEDgC/A9YDp2OSwoOYtv/O1vMgTE+hEZimn1eBvhx9FaDRQEXEew4dgqeegrlz\nYfhwmDPHFHpbKNbRQOO9AngPeBrYABwB3gGeAE4GVgJTMcXeKdb2W6z1W4DDwM2oCUhExAgGYdUq\n+MtfIDMz6X9O8wGIiLiE5gMQEUl1kQh89ZXdUSgBiIi0mUgEXnrJNO/MmWN3NBoLSESkTaxbB/fe\nC3v2wP33w6WX2h2REoCISFJFIjBxImzdCr/8JVx1FbRvb3dUgIrAIiLJt349DBsGPl9S/0ysRWAl\nABERl1AvIBERO5SWwrx5dkcREyUAEZF4BINw441wwQWQlgY1NXZH1GJKACIirREKwc03myEbevaE\nsjK4556UKfC2hHoBiYi0Rn4+dOxoJmPp0sXuaFpFRWAREZdQEVhEJJH27YMjR+yOIimUAEREGlM7\nGUvfvlBcbHc0SaEEICIS7cABWLjQnPjfeQdefx3OO8/uqJJCRWARkVrbtsHFF8M558Dq1TB0qN0R\nJZWKwCIitaqr4b334Nxz7Y6kVTQUhIiIR6kXkIhIcyIReP55ePZZuyOxnRKAiHhDJGLa9UeMMMMy\nd+pkd0S2UxFYRNwvEDCTsXzxhTn5X345tNPvX9UARMTdIhEzCcuECXD11Y4aqydWKgKLiHiUisAi\n4l179tgdgaMoAYiI85WVwTXXwMiRjhqP325KACLiXMEgTJ0K558PgwaZMXtc3MafaEoAIuJMixeb\nyVh69DBXADNnwskn2x2Vo6gILCLO9MEH5oR/2ml2R5Iy1AtIRFJaOBymoKCI/PxCKir207t3B3Jz\nsxg1aiTp6el2h+doSgAikrKmTJlFcfEJhEJZVFePBNKBMD5fEX5/IZmZB1m58v66D+zbB48+aiZd\n79nTrrAdw45uoJ2BfGArsAU4DzgVWAOUAq9Y29SaAZQB24BxCfj7IuIA4XCY4uITCAZnUl09FnPy\nB0inunosweBMiouPJxwO15+Mpbwc0lL1t6qzJSIBLAT+DgwEhmBO7HdjEsCZwFrrNcAg4ErrOQdY\nlKAYRCTFFRQUEQplNbvN7srhfDTtTnPif/ttM4TD00+D3982QXpMvCffTsD3gKes14eBvcBEYJm1\nbhkw2VqeBKwAqoEgUA6MiDMGEXGA/PxCq9mnad0O+znw0jp46SVYudJ07ZSkiTcBnAF8BvwWeAdY\nApwEdAWqrG2qrNcAPYDKqM9XAkrtIh5QUbGfumafxpUzhLv6XwZnn902QXlcvKOBHgecA9wKFAOP\nUtfcUytiPZrS6Ht5eXnfLmdnZ5OdnR1HmCJit969OwBhIJ00jnAK+9hbrzwIELa2k5YIBAIEAoFW\nfz7eyko3oBBzJQBwIabI+2/AGGAX0B1YBwygLjk8YD2vBmYDbzX4XvUCEnGZ1atfY+KECOMPf80c\n7uNv/Af38Ot62/h8a1m1Ko2cnLE2RelssfYCivcKYBewE1PsLQUuAkqsx/XAg9bz89b2q4DlwCOY\npp9+wPo4YxCRVBeJkH1gH8XtbiVCF2Yyjxe55KjN/P5CRo2aZkOA3pSICWF+DjwDHA/sAG4A2gMr\ngamYYu8Ua9st1votmILxzTTfPCQiTheJQE4OJ+7cydpzRvLYJ0Oo/DgdqmtrAmF8vkLrPoBDuhms\nDaVq51o1AYm4yTvvwNCh0L697gROIt0JLCLiUZoQRkTsUVICeXmmyUccQQlAROJTVgbXXgtjx5rR\nOZUAHEMJQERaJ3oyloEDzZg9//Vf0E6nFadIRC8gEfGiF1+E7t2htBS+8x27o5FWUBFYxOPUK8c9\n1AtIRFqsRePzP/4L6NRJTTsOoF5AItIixxqf//Pgzxm55k0i/fqZoZnFdZQARDyqqfH50/mGO5lP\nGf34173teGP+AsjMtCFCSTYlABGPamx8/gFspZy+ZFLMGNZxXeSvPP3PCpsilGRTLyARj2psfP4y\n+pHDajYxtMF24ka6AhDxqLrx+evUcFy9k7/G53c3JQARrzlyBP70J2497QA+X1Gzm/p8heTmNj+P\nrziXEoCIV0Qi8Ne/mukWH3mEAaMvwO8vbPYjZnz+5ufxFedSDUDE7SIRePlluO8+qK6GefPgkks4\nMS2NzMx/AnOt+wCy0Pj83qIbwUTcLhKB666DSZPg8suPuqFLdwK7h+4EFhHxKN0JLOJlu3fbHYE4\niBKAiBuUlEBuLowaBTU1dkcjDqEEIOJktZOxjBkD550HGzZA+/Z2RyUOoQQg4lSPPQZZWdC/v5mM\n5c47QUVbiYGKwCJOFQzCKafAqafaHYmkCPUCEhHxKPUCEnGTL7+Ee++FCo3IKYmnBCCSivbtgzlz\noF8/qKoCn8/uiMSFlABEUkk4DA89ZE78paVQVARLlpjJ10USTGMBiaSSTz6B4mJ47TUYPNjuaMTl\nVAQWEXEJu4rA7YF3gRes16cCa4BS4BWgc9S2M4AyYBswLkF/X8RZamrgiy/sjkI8LlEJ4HZgC1D7\ns/1uTAI4E1hrvQYYBFxpPecAixIYg0jqO3IE/vxn+O534ZFH7I5GPC4RJ9+ewL8DS6m79JgILLOW\nlwGTreVJwAqgGggC5cCIBMQgktoiEVi1CoYNM0XehQtNLx8RGyWiCLwAuBM4JWpdV6DKWq6yXgP0\nAKLnoKsE/AmIQSR1RSLw/e/D55+bk/6ECZCWquU38ZJ4E8AlwKeY9v/sJraJUNc01NT7R8nLy/t2\nOTs7m+zspr5eJMWlpZlf/IMHHzUZi0g8AoEAgUCg1Z+P92fIr4DrgMPAiZirgOeATExC2AV0B9YB\nA6irBTxgPa8GZgNvNfhe9QISEYmRnWMBjQb+G5gAzAc+Bx7EnPQ7W8+DgOWYdn8/8CrQl6OvApQA\nJGW0eMrEd981Bd5589TEI7aINQEk+kaw2rP2A8BKYCqm2DvFWr/FWr8Fc9VwM803D4nYasqUWRQX\nn2BNmv4LaidNf/rpIvz+BWRmHmRl3lUweza8+SbMmGHa/JUAxAFS9V+prgDEduFwmMGDFxAMzmz0\n/T6U89BJuUzuECJt+nS45RaNxy+2svsKQMQ1CgqKCIWymnz/Il5l0/6z6bhsDhdfPqENIxNJDHVJ\nEGlCfn4h1dUjm3z/cX5G3pFF/OmlTW0YlUjiKAGINKGiYj+Qzr+wm3Y0NdF6urWdiPMoAYg0YUBX\nmMNdbKc/w3i3ia3C9O7doU3jEkkUJQCRhqzJWP5n1f/hT9vEcN7mbc5tdFOfr5Dc3KbrBCKpTAlA\nJFpJCfTtC9u3U/PGG9x/+oV8REaTm/v9hYwa1XSdQCSVqReQSLT+/eH112HgQDoAmZn5wFzrPoAs\nau8D8PkK8fsLycw8VP9mMBEH0X0AIsfQ4juBRWxm51AQiaQEIMlTUwPLl0P79vDDH9odjUjC2DUj\nmEjqq52M5ayz4Ikn4PTT7Y5IxFaqAYj7RSLwwgswaxb4fLBgAYwbp/F6xPNS9X+AmoAkcSIRmDoV\nJk2CiRN14hfXUg1ARMSjVAMQb6uqOvY2IgIoAYhbbNxomnfGjjXFXhE5JiUAcbYtW+CKK2D8eLjo\nInj7bc27K9JC+p8izvXoo5CdDZmZUF4Ot90GJ55od1QijqEisDjXzp3QqROccordkYikBPUCEhHx\nKPUCEnf57DOYPh0++sjuSERcRwlAUtOePXDffTBgAHz9tdr2RZJACUBSy1dfwbx50K8ffPyx6dWz\naBF07Wp3ZCKuo7GAJG4JHS75s89g61Z4800488zkBCwigIrAEqcpU2ZRXHyCNWHKSOomTCmyJkw5\nyMqV99sdpognxFoE1hWAtFo4HKa4+ASCwZkN3kmnunosweBYYC7hcPjoK4Hqati7F7p0aatwRaQB\n1QCk1QoKigiFmp8QPRTKoqCgqG5FTQ384Q8wcCA89liSIxSR5igBSKvl5xdazT5Nq67OIj+/sP5k\nLIsXw5IlkJfXNoGKSKPUBCStVlGxH9Pm35x0Kj4Kw+jRsH8/PPII/OAHGpNfJAXEewXQC1gHlADv\nA7dZ608F1gClwCtA56jPzADKgG3AuDj/vtiod+8OQPgYW4XpfXq6+cVfXAw5OTr5i6SIeBNANfAL\nYDAwErgFGAjcjUkAZwJrrdcAg4ArreccYFECYhCb5OZm4fMVNbuNz1dIbm6WuaFLJ36RlBLvyXcX\nsNFa/hrYCviBicAya/0yYLK1PAlYgUkcQaAcGBFnDGKTUaNG4vcXfvs6k/XM506grguv31/IqFHN\n1wlExB6J/PWdAQwD3gK6ArVTM1VZrwF6AJVRn6nEJAxxoPT0dDIzDzK+xy28kHY+z3IZO+hDGt/g\n860lI2MumZmHYr8ZTETaRKKKwB2BZ4Hbga8avBch+ifh0Rp9Ly+qh0h2djbZ2dlxBShJsG0bK9O2\nE6l5nW0/uZx5B8fxwceV3Nh7oXUn8DSd/EWSKBAIEAgEWv35RDTK+oAXgZeAR61124BsTBNRd0yh\neAB1tYAHrOfVwGzMVUM03QnsBEuWwBdfwK23wkkn2R2NiOe19XwAaZg2/s8xxeBa8611D2JO+p2t\n50HAcky7vx94FejL0VcBSgAiIjFq6wRwIVAAbKLuJD4DWA+sBHpjir1TgD3W+/cANwKHMU1GLzfy\nvUoAqWTXLjjtNGjf3u5IRKQZmhFMEmf3bpg/H5YuhTVrYPhwuyMSkWZoRjCJ3549MGsW9O9vJmPZ\nvFknfxEX0lAQUt/778OYMTBhAmzYAGecYXdEIpIkagKS+mpq4MMPoW9fuyMRkRipBiAi4lGqAcix\nVVfDk0/C739vdyQiYiMlAC+pqYFnnoFBg2D5cs25K+JxKgJ7QSQCzz1nevZ07gxPPGEKvSLiaUoA\nXvHKK/Dww5qMRUS+lapnAhWBRURipCKw14VCdkcgIg6hBOAWtdMtjh9vJmAXETkGJQCne+89mDQJ\nLr0UJk82d++202EVkWPTmcLJaou6Y8ZAWRn87Gdw/PF2RyUiDqEisJOFQtCpE3TsaHckIpICNBSE\niIhHqReQ2+zaBbffDsGg3ZGIiMsoAaSq3bth+nQzbEO7dmrmEZGEUwJINXv3wuzZ9SdjWbAAunSx\nOzIRcRkNBZFqvvwSKis1GYuIJJ2KwCIiLqEisFMcPAhVVXZHISIepgTQ1g4fNpOxnHkmLFlidzQi\n4mGqAbSVmhr44x8hLw9694YVK+D88+2OSkQ8TAmgLUQicOGFZhz+xx+HsWPtjkhEREXgNlNeDn36\naDIWEUkaDQUhIuJR6gVkp3/8A+64wzT5iIikOCWARNiwwUzEcu21MGSIEoCIOIJdCSAH2AaUAXc1\ntsGPfzyP1atfIxwOt2lgMdm82UzCMnkyTJwIpaVw442akEVEHMGOGkB7YDtwERACioGrga1R20Tg\nG3y+Ivz+QjIzD7Jy5f02hHoMv/2tGbvnpz+FDh3sjkZEPM4JReAsYDbmKgDgbuv5gahtIlDXjJKR\nMZeSkmmkp6e3TYQiIg7khCKwH9gZ9brSWtekUCiLgoKipAbVrMpKcyOXiIiL2JEAYq6QVldnkZ9f\nmIxYmlc7GcuQIbBpU9v/fRGRJLLjTuAQ0CvqdS/MVUADeVHL2VRU7E9qUPV8/jnMnw9Ll8J118GW\nLdCtW9v9fRGRFggEAgQCgVZ/3o4awHGYIvD3gY+B9TRaBI6+UAgzdeoCli6dmfzoNm+GMWPgiitg\n5kzo2TP5f1NEJAFirQHYcQVwGLgVeBnTI+hJ6p/8j+LzFZKbm9UGoWGmYNywATIy2ubviYjYJGWH\nglAvIBGR2DjhCqCFwvh8hdZ9AIcSe/I/dMi073foADfckLjvFRFxkJS9ZXXq1AWsWpVGScm0xN0E\ndvgwPPWUmYzlhRdM7x4REY9K2SaghI4GGomYCVjy8sDvh7lz4YILEvf9IiIpwEVNQAlWWAiLF2sy\nFhERizeuAEREPMAJQ0EkTyQCwaDdUYiIOIJ7EsCbb5rmncsugyNH7I5GRCTlOT8B1E7Gcs018KMf\nwfr1Go9fRKQFnH2mfPDB+pOx3HADHOeduraISDycXQTetQs6ddJkLCIiOGNCmJZQLyARkRi5rxdQ\nZSXcdBN88IHdkYiIuErqJoBdu+COO2DoUNPM07mz3RGJiLhK6lZMBw82k7GUlGgyFhGRJEjdGkBF\nBfTqdewtRUQEUBFYRMSz3FcEFhGRpFACEBHxKCUAERGPUgIQEfEoJQAREY9SAhAR8SglABERj1IC\nEBHxKCUAERGPUgIQEfEoJQAREY9SAhAR8ah4EsBDwFbgPeA5oFPUezOAMmAbMC5q/XBgs/Xewjj+\ntoiIxCmeBPAKMBgYCpRiTvoAg4ArreccYBF1o9P9BpgK9LMeOXH8fccKBAJ2h5A0bt430P45ndv3\nL1bxJIA1wBFr+S2gp7U8CVgBVANBoBw4D+gOnAyst7Z7Gpgcx993LDf/I3TzvoH2z+ncvn+xSlQN\n4Ebg79ZyD6Ay6r1KwN/I+pC1XkREbHCsKSHXAI3Nx3gP8IK1PBM4BCxPYFwiIpLi/hN4Ezgxat3d\n1qPWakwTUDdM0bjW1cDiJr63HIjooYceeugR06OcNpIDlABdGqwfBGwEjgfOAHZQVwR+C5MM0jBN\nRp4sAouIOF0Z8BHwrvVYFPXePZhMtA34QdT62m6g5cD/tk2YIiIiIiKSsnIwVw1lwF02x5IMQWAT\n5oppffObOsJTQBXmqq7WqZjOA6WYe0U62xBXojS2f3mY3my1V75ObcbsBazDNOO+D9xmrXfL8Wtq\n//Jwx/E7EdOkvhHYAvzaWu/Y49ce0zSUAfgwOzbQzoCS4EPMAXKL7wHDqH+CnA9Mt5bvAh5o66AS\nqLH9mw1MsyechOoGnG0tdwS2Y/6/ueX4NbV/bjl+AOnW83FAEXAhMR6/VBoLaAQmAQQxN5H9EXNT\nmdukHXsTx3gD+LLBuonAMmt5Gc6+2a+x/QN3HMNdmB9ZAF9jeuj5cc/xa2r/wB3HDyBsPR+P+QH9\nJTEev1RKAH5gZ9Tr2hvI3CQCvApsAH5icyzJ0hXTbIL13NXGWJLl55gxsJ7EQZfYzcjAXOm8hTuP\nXwZm/4qs1245fu0wSa6KuuaumI5fKiWAiN0BtIELMP8QxwO3YJoY3Ky2b7Kb/AbTvfls4BPgYXvD\niVtH4FngduCrBu+54fh1BPIx+/c17jp+RzD70RMYBYxp8P4xj18qJYAQpnBTqxf1h45wg0+s58+A\nv2Cavdymirq7x7sDn9oYSzJ8St1/rKU4+xj6MCf/3wPPW+vcdPxq9+8P1O2fm45frb3A3zDd7GM6\nfqmUADZgRgjNwLRpXQmssjOgBEvHDIYHcBJmmOzNTW/uWKuA663l66n7j+cW3aOWL8W5xzAN0wSy\nBXg0ar1bjl9T++eW49eFuuarDsDFmF5Njj5+4zHV+nLqhpd2izMw7XUbMd3S3LB/K4CPMWNB7QRu\nwPRyehUHdkNrRMP9uxEziu0mTBvy8zi3jfxCTBPCRup3iXTL8Wts/8bjnuN3FvAOZv82AXda691y\n/EREREREREREREREREREREREREREREREREREBOD/AY6nw+dT8vwkAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10c1f9350>"
       ]
      }
     ],
     "prompt_number": 54
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example from the textbook: quadratic (second order polynomial) interpolation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline\n",
      "\n",
      "def evalPoly(a,xData,x):\n",
      "    n = len(xData)-1  # order of polynomial\n",
      "    p = a[n]\n",
      "    for i in range(1,n+1):\n",
      "        #print \"a[\", (n-i),\"]=\", a[n-i]\n",
      "        p = a[n-i] + (x-xData[n-i])*p\n",
      "    return p\n",
      "\n",
      "def coefficients(xData,yData):\n",
      "    n = len(xData)\n",
      "    a = np.zeros((n, n));\n",
      "    for i in range(0, n):\n",
      "        a[i, 0] = yData[i];\n",
      "    for j in range(1, n):\n",
      "        for k in range(0, n-j):\n",
      "            a[k, j] = (a[k+1, j-1] - a[k, j-1])/(xData[k+j]-xData[k])\n",
      "    return a[0, :] # return the zeroth row\n",
      "\n",
      "xAll = [0, 10, 15, 20, 22.5, 30]\n",
      "yAll = [0, 227.04, 362.78, 517.35, 602.97, 901.67]\n",
      "plt.plot(xAll,yAll,'bo',markersize=10)\n",
      "\n",
      "xData = xAll[1:4]\n",
      "print \"xData=\",xData\n",
      "yData = yAll[1:4]\n",
      "print \"yData=\",yData\n",
      "a = coefficients(xData, yData)\n",
      "print \"a=\",a\n",
      "y = evalPoly(a, xData, 10)\n",
      "print \"f(10)=\", y\n",
      "y = evalPoly(a, xData, 15)\n",
      "print \"f(15)=\", y\n",
      "y = evalPoly(a, xData, 16)\n",
      "print \"f(16)=\", y\n",
      "\n",
      "xs=np.linspace(0, 30, 50)\n",
      "ys=np.zeros(50)\n",
      "i=0\n",
      "for x in xs:\n",
      "    ys[i] = evalPoly(a,xData,x)\n",
      "    i=i+1\n",
      "plt.plot(xs,ys,'r--')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "xData= [10, 15, 20]\n",
        "yData= [227.04, 362.78, 517.35]\n",
        "a= [ 227.04     27.148     0.3766]\n",
        "f(10)= 227.04\n",
        "f(15)= 362.78\n",
        "f(16)= 392.1876\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 56,
       "text": [
        "[<matplotlib.lines.Line2D at 0x10bd38c90>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGddJREFUeJzt3Xt4VeWd6PEvYkCC9TZW0GBOGMQLTq03IIwa441Sx/HS\n4eC9jJeOHltrtbbqYR6hDq23Wi/njHXa6gitQmmsaD3qUak5URsgKCByUaATkShBFATcCIHs88da\nmBCSQHaSvfZa+/t5njx77bXX3vu3ngW/vPm973pfkCRJkiRJkiRJkiRJkiRJUg57DKgHFjTbdwDw\nMvAe8BKwX7PXbgOWAkuAkc32nxB+xlLgwW6MV5LUSacAx7Fj4r8H+HG4fQtwV7g9BJgHFAAlwDKg\nR/jabGBYuP08MKrbIpYkdVoJOyb+JUC/cLt/+ByC1v4tzY57ESgFDgYWN9t/EfBIdwQqSdq1PTJ4\nTz+C8g/h4/ZfAocAK5sdtxIoamV/XbhfkhSBTBJ/c+nwR5IUE3tm8J56ghLPKoIyzupwfx1waLPj\nBhC09OvC7eb761r74EGDBqWXL1+eQUiSlNeWA4ft7sGZtPifBcaG22OB6c32XwT0AgYCgwk6dVcB\n64HhBJ29lzd7zw6WL19OOp1O7M/48eMjj8Fz8/w8v+T9AIM6ksR31eKfApwKHAh8ANxOMIpnGnAV\nUAuMCY9dFO5fBGwFrqOpDHQd8DjQh2BUz4sdCVKS1HV2lfgvbmP/mW3s/1n409KbwNd2NyhJUvfp\nbOeuOqC8vDzqELpNks8NPL+4S/r5dVSPXR+SVemwXiVJ2k09evSADuRzW/ySlGdM/JKUZ0z8kpRn\nTPySlGdM/JKUZ0z8kpRnTPySlGdM/JKUZ0z8kpRnTPySlGdM/JKUZzJZiEWSFLFUKkVV1UwqKqo7\n/F4naZOkmBkz5nZqanpTVzeChoZSoC90IJ/b4pekGEmlUtTU9Ka2dlzGn2GNX5JipKpqJnV1Izr1\nGSZ+SYqRiorqsLyTORO/JMXIihWbgMJOfYaJX5JipLi4D5Dq1GeY+CUpRkaPHkFBwcxOfYaJX5Ji\npKyslKKiYOz+AD7gYp7s8GeY+CUpRgoLCxk6dDNj+l/NLI7jYN7v8GeY+CUpTtJpppUdxNRtz1I3\n8cd8dlVjhz/CO3clKU6mTYOJE+Hpp2HQIAB69OgBHcjnJn5JipOtW2HzZujb98tdJn5JyjMdTfzW\n+CUpz5j4JSkXbdoEN90E9fVd/tEmfknKNe+/DyefDKtWwd57d/nHm/glKZfMmAHDh8Nll8ETT+zQ\nidtVnI9fknLFfffBz38OU6bAaad129eY+CUpV+y7L8yaBcXF3fo1DueUpJhzOKckqV0mfknKtm3b\nYPnyyL6+M4n/NmAhsAB4EugNHAC8DLwHvATs1+L4pcASYGQnvleS4uvjj2HUKBg/PrIQMk38JcB3\ngOOBrwE9gYuAWwkS/+HAjPA5wBDgwvBxFPBwJ75bkuKppgZOPBGGDoVJkyILI9Pkux5oIFj4cc/w\n8UPgXGD72UwCzg+3zwOmhO+pBZYBwzL8bkmKl3QaHnkE/uEf4IEH4Gc/g549Iwsn0+GcnwL3ASuA\nTcD/JWjp9wO2319cHz4HOARovlbYSqAow++WpHhZsCBI/G+8AYMHRx1Nxol/EPADgpLPZ8AfgMta\nHJMOf9rS6msTJkz4cru8vJzy8vIMQ5SkHHHMMfDmm13Wyq+srKSysjLj92c6jv9C4Czg6vD55UAp\ncDpwGrAKOBh4FTiSplr/XeHji8B4YFaLz3UcvyR1ULbG8S8hSPR9wi87E1gE/AkYGx4zFpgebj9L\n0PnbCxgIDAZmZ/jdkpS7YtB4zTTxzwcmA3OAt8N9vyJo0Z9FMJzzdJpa+IuAaeHjC8B1tF8GkqT4\n+egjOPNMmDMn6kja5ZQNktQVKivh0kvhmmtg3LisjtrpaKnHSdokqTMaG+Huu+Ghh2DyZDjrrKgj\n2iUTvyR1xtVXw7vvBjdnDRgQdTS7xVKPJHXG4sVw2GFQUBBZCB0t9Zj4JSnmnJZZktQuE78k7Y75\n84NpFxLAxC9J7Umn4de/Dsbn77NP1NF0CUf1SFJbNm6Ea68NWvuvvQZHHhl1RF3CFr8ktebdd4N5\n83v1ChZAT0jSB0f1SFLrPvwQ/vxnuKzlxMO5x+GckpRnHM4pSWqXiV9Sfkun4ZlnYNu2qCPJGkf1\nSMpf69cHs2kuWgSlpdCv367fkwC2+CXlpzlz4PjjYf/9YebMvEn6YOKXlG/Sabj/fjj7bLjrLnj4\nYejTJ+qosspSj6T80tgItbXB2PyBA6OOJhIO55SkmHM4pySpXSZ+Scm1ciXU10cdRc4x8UtKpunT\n4YQTgkXQtQM7dyUly6ZNcPPN8PzzQfIfMSLqiHKOLX5JybFoEQwfDmvWwNy5Jv02OKpHUnLcfTcc\neCBceSX0yLX01n2cnVOS8ozDOSVJ7TLxS4qfhoZghSxlxMQvKV6WL4eTTw7q+cqIiV9SPKTTMGlS\nMH3ypZfCo49GHVFsOY5fUu5buxauvRYWLoQZM+CYY6KOKNZs8UvKffPmBfPl19SY9LuAwzklKeYc\nzilJapeJX1LuSKdh9uyoo0g8O3clZUUqlaKqaiYVFdWsWLGJ4uI+jB49grKyUgoLC2HVKrj6avj4\nY3j9dSgoiDrkxOpMjX8/4DfA0UAauAJYCvwe+G9ALTAGWBcefxtwJbAN+D7wUiufaY1fSqAxY26n\npqY3dXUjaGgoBQqBFAUFMykqquZ/HPI2P/5rFVx1FYwfb9LvoGzO1TMJ+H/AYwR/OfQFxgFrgHuA\nW4D9gVuBIcCTwFCgCHgFOBxobPGZJn4pYVKpFEcffT+1teN2eq0vG7mfG/nGnk9x0It/YK8zzogg\nwvjLVufuvsApBEkfYCvwGXAuwS8Ewsfzw+3zgClAA8FfAsuAYRl+t6QYqaqaSV1d69Mj9+VzNtGH\n43icyoZcG2SYXJkm/oHAx8B/Am8BvyZo8fcDtq9zVh8+BzgEWNns/SsJWv6SEq6iojos7+xsNf24\ngYf4dOuZVFRUZzmy/JVp4t8TOB54OHz8nKCk01w6/GmLNR0pD6xYsYmgpt+ewvA4ZUOmo3pWhj81\n4fMKgs7bVUD/8PFgYHX4eh1waLP3Dwj37WTChAlfbpeXl1NeXp5hiJJyQXFxH3qwkYt5hqlcRCM9\nWzkqRXFxn6zHFleVlZVUdmIt4c4U1aqAq4H3gAk0/Ur/BLib4C+A/dixc3cYTZ27h7Fzq9/OXSlh\nKidPIf3P91CQ7sv5TOcTDtzpmIKCGTz7bA9GjTo9ggjjr6Odu50Zx3898ATQC1hOMJyzJzANuIqm\n4ZwAi8L9iwg6gq/DUo+UbOk0/Pa3nPrDH3LPvsfxP9e90EZrH4qKqikruynLAeavXOtGt8UvJcGG\nDfDtbwdz50+ezJif/ZGaml7hOP4RNI3jr6aoqJqhQ7cwbdodUUcdW665Kyl6jY3wq1/BFVdA797A\nbty5q4yZ+CUpzzg7pySpXSZ+SZnbsAFuvBE+/DDqSNQBJn5JmXn11WA1rPXroW/fqKNRBzgts6SO\n+fxzuPVWePrpoAP37LOjjkgdZOKXtPs2b4YTToDhw2HBAth//6gjUgYc1SOpY959F444Iuoo1IzD\nOSUpz2RzygZJMbbLG6oaGlwJK6Fs8Ut5aFdLIX675D1+svIvUFEBX/961OFqF2zxS2pXKpWipqZ3\nK0shFlLQMJwba6cz5oOn2fzkb+ht0k8kx/FLeaatpRBPpZL5fJ19+Yyv95jEq/vsPH2yksHEL+WZ\n1pZCLGALP2UcP+AB/plJrN76DZdCTDBLPVKeaW0pxAZ6cTKv01QmdinEJLPFL+WZYInDVCuvNO8b\ndCnEJDPxS3nm2uI96LtnVbvHFBRUM3r0zv0ASgYTv5QvVq2C0aM5/onHOb7fy+0eGiyFWNruMYov\nE7+UdOk0PPZYMJPmEUewx/z59P/7vpSUTKSgYAZNZZ8UBQUzKCmZyNChW1wVK8G8gUtKso0b4YIL\nYO1aePTRHW7GcinE5HCuHklN0ml46ik4/3zY00F8SWXil6Q845q7Ur6y0aTdZOKXkmDGjKB+X18f\ndSSKAYt+Upx98gncfDP8+c/w7/8O/fpFHZFiwBa/FEfpNEyZAn/3d7DPPvDOO3DOOVFHpZiwxS/F\n0cqV8OCDMH16sP6t1AGO6pHiKp2GHrn2X1hRcFSPlC9M+sqQiV/KZevWwS9/6VBNdSkTv5SL0mmY\nOhWOPhrmzYOtW6OOSAli566Ua5Yuhe9+NxiTX1EBI5weWV3LFr+USyorg0Q/ahS8+aZJX90i13qH\nHNWj/PbFF/Dxx3DooVFHohhxkjZJyjMO55TiYOtW+Otfo45Ceaqzib8nMBf4U/j8AOBl4D3gJWC/\nZsfeBiwFlgAjO/m9UnzNnAnDhsFPfhJ1JMpTnU38NwCLgO31mVsJEv/hwIzwOcAQ4MLwcRTwcBd8\ntxQvn3wC//Iv8K1vBROrPf541BEpT3Um+Q4AzgZ+Q1Nt6VxgUrg9CTg/3D4PmAI0ALXAMmBYJ75b\nipeKChgyBPbaCxYvhksu8c5bRaYz4/jvB34E7NNsXz9g+4Tg9eFzgEOAmc2OWwkUdeK7pXjp3x9e\nfBGOOy7qSKSME/85wGqC+n55G8ekaSoBtfW6lB9OPjnqCKQvZZr4/56grHM2sBdBq/+3BK38/sAq\n4GCCXw4AdUDzgckDwn07mTBhwpfb5eXllJeXZxiiFIHGxmDETq9eUUeiBKusrKSysjLj93dFkfFU\n4GbgH4F7gE+Auwk6dvcLH4cATxLU9YuAV4DD2LnV7zh+xdfs2fC978GVV8K110YdjfJIR8fxd9Vc\nPduz9V3ANOAqgk7cMeH+ReH+RcBW4Dos9SgpVq+G226DF16AO++Eyy+POiKpXbk2rMAWv3JGKpWi\nqmomFRXVrFixieLiPowePYKyslIKCwuDGTQfeggmToSxY+H224NlEKUsc8oGqQuMGXM7NTW9qasb\nQUNDKVAIpCgomElRUTVDh25m2rQ7gqT/T/8ERx0VdcjKYyZ+qZNSqRRHH30/tbXj2jympGQiCxfe\nFLT8pYg5V4/USVVVM6mr23E65D3YtsPzuroRVFXNRIojE7/UQkVFdVjeAUhzIVNZwpH0Y9WXxzQ0\njKCiojqaAKVOcgUuqYUVKzYBhRzLXB7kBr7CBq7iUerp3+yowvA4KX5s8UstDPlqI49wFS/wTX7H\nZZzIHF6jrMVRKYqL+0QSn9RZtvilFs4f+TXmTl3KUY2LWcf+rR5TUFDN6NEui6h4ssUvtTDsv5/H\nQ8XHtpn0AYqKqikrK23zdSmX2eJXftu2DXr23GFXYWEhQ4duBiaG4/hH0DSOvzocx7/FoZyKLcfx\nKz+tWhXcabtxIzz5ZKuH7PLOXSlHeAOX1J5Nm+D+++EXv4ArroBx42C//Xb9PimHRTVJm5T7Kirg\nhz8M1rudNQsGDYo6IikSJn7lj02b4IknXBRFec9SjyTFnHP1SGvXBlMmS2qViV/J8cUXcO+9cPjh\nMH9+1NFIOcvEr/hrbAxq90ceCW+8Aa+/DsceG3VUUs6yc1fxtmIFXHBBcBPW5MlQ1nJOHUkt2bmr\neNuyBZ5/Hs47D3rk2j9nKTu8gUuS8oyjepRMn34Ks2dHHYWUCCZ+5bZUCu68E444Ap57LupopEQw\n8Ss3NTTAf/wHDB4Mc+cGo3XuuCPqqKREcFSPctPFF8O6dfDMM3DiiVFHIyWKnbvKTevWOWumtJsc\n1SNJecZRPYqPt9+GSy4JWveSssbEr+xbtgwuvRRGjoThw6FPn6gjkvKKiV/Zs3IlfOc7UFoazKuz\ndCnccAP07h11ZFJecVSPMtbhNWnr6uCgg+C99+CAA7IfsCTAzl1laMyY26mp6U1d3QgaGkqBQiBF\nQcFMioqqGTp0M9OmOe5eygbX3FW3S6VS1NT0prZ2XItXCmloOJ11tcexdNu9pFKp1lv+kiJljV8d\nVlU1k7q6ETvt34fPuJ2fsJTBHPNhcJyk3GPiV4dVVFSH5Z3AV1jPOCayjMMYyH8xnFlM3vavVFRU\nRxilpLZY6lGHrVixiaCmD3uzgSUcyQzO4CTeYCmHtzhOUq4x8avDiov7ACmgkI18heN5i3r6tzgq\nFR4nKddkWuo5FHgVWAi8A3w/3H8A8DLwHvAS0HyylduApcASYGSG36scMHr0CAoKmur3Oyd9KCio\nZvTonfsBJEUv08TfANwIHA2UAt8FjgJuJUj8hwMzwucAQ4ALw8dRwMOd+G5l2/r18NOfwnXXAVBW\nVkpRUfv1+6KiasrKSts9RlI0Mk2+q4B54fZGYDFQBJwLTAr3TwLOD7fPA6YQ/MKoBZYBwzL8bmXL\nunXBHPiDBsHixfD94A+7wsJChg7dTEnJRAoKZhCUfSAYxz+DkpKJDB26xaGcUo7qihp/CXAcMAvo\nB9SH++vD5wCHAM3H9q0k+EWhXHXnnXDffXDOOfD668EKWM1Mm3ZHszt3729x5+5NJn0ph3U28e8N\nPAXcAGxo8Vo6/GlLq69NmDDhy+3y8nLKy8s7FaAyNGhQsMbt3/5tm4cUFhYyatTpjBp1ehYDk1RZ\nWUllZWXG7+/MlA0FwHPAC8AD4b4lQDlBKehggg7gI2mq9d8VPr4IjCf4K6E5p2yQpA7K1nz8PYBH\ngUU0JX2AZ4Gx4fZYYHqz/RcBvYCBwGBgdobfra6yfDnce2/UUUjKskwT/0nAZcBpwNzwZxRBi/4s\nguGcp9PUwl8ETAsfXwCuo/0ykLrTwoVw2WXBXPjr18PWrVFHJCmLnJ0zn7z1Fvzbv0F1NfzgB8Hw\nzH32iToqSZ3k7Jxq25tvwmmnwRNPgKNupLxli1+SYs7F1vPdli0wdSo0NkYdiaQcZeJPivXr4ec/\nD8bdP/YYfPpp1BFJylEm/rj76CO49VYYODCo4f/pT/DSS3DggVFHJilH2bkbd6++Cp9/DnPmBMlf\nknbBzl1Jijk7d5No61aoqIAvvog6EkkJYOLPZRs2wIMPwuDBweOqVVFHJCkBTPy5qK4ObrklqNm/\n8UYwPPO116CkJOrIJCWAnbu5aNEi2Lx5l9MiS1Im7NyVpJizczcu1q2DX/wC1qyJOhJJecbEn23v\nvgvf+15QwnnzTUildv0eSepCJv5smTcPzj4byspg//3hnXeCWTKLi6OOTFKesXM3W/bYA0aPhj/+\nEfbaK+poJOUxO3clKebs3I3Ktm3BBGnf+AYsWxZ1NJLUJks9nbVmDTz6KDzyCHz1q3D99TBgQNRR\nSVKbbPF3xtSpwXQKS5bA738f3HB1+eXW8CXlNGv8nbFmDfToAX/zN1FHIimPWePvDkuXQmu/kA48\n0KQvKXZM/G354otgnP0pp8Cpp0J9fdQRSVKXMPG3tHw5/OhHcOihMGkS3HgjvP8+9O8fdWSS1CUc\n1dPSK68Ej9XVcNhh0cYiSd3Azl1Jijk7d3fl88+DEs6YMdDYGHU0kpR1+ZH402mYNQuuuSao3f/h\nD3Dxxa2P1JGkhMuPGv/YsfCXv8CVV8KCBVBUFHVEkhSZ/Kjxr14djLnfIz/+wJGUX/K3xr94MTz7\nbOuvHXSQSV+SQvHOhmvXBpOjlZbCGWfAwoVRRyRJOS+epZ7GRrj0UnjhBTjrLLjiChg5EvbMjy4L\nSWquo6WeeCZ+gOeeg5NOCpYxlKQ8lqzE/8EH0NAQLEwuSWpVrnfujgKWAEuBW1o74PrL/5UFN97M\ntlNPhWOPhddey2qAkpR02Uz8PYH/TZD8hwAXA0e1PGji7/4Xyx6czXWLB3BJ+TXBGPyEqKysjDqE\nbpPkcwPPL+6Sfn4dlc3EPwxYBtQCDcBU4LyWB5VQy7fSVfzq4yeofquQVCqVxRC7V5L/8SX53MDz\ni7ukn19HZTPxFwEfNHu+Mty3g3U0ddbW1Y2gqmpm90cmSXkkm4m/w7fkNjSMoKKiujtikaS8lc1R\nPaXABIIaP8BtQCNwd9Mhg9KwPIshSVIiLAdycgGRPQmCKwF6AfNopXNXkpQs3wTeJejkvS3iWCRJ\nkiRlyy5v7Iq5WuBtYC4wO9pQusRjQD2woNm+A4CXgfeAl4D9Ioirq7R2fhMIRqLNDX9G7fy2WDgU\neBVYCLwDfD/cn5Tr19b5TSAZ128vYBZBqXwRcGe4P3bXrydB6acEKCCZtf//IrgwSXEKcBw7JsZ7\ngB+H27cAd2U7qC7U2vmNB26KJpwu1R84Ntzem6D0ehTJuX5tnV9Srh9AYfi4JzATOJkOXr9cmJZ5\nt27sSoBcmxepM14D1rbYdy4wKdyeBJyf1Yi6VmvnB8m4hqsIGlcAG4HFBPfTJOX6tXV+kIzrB7D9\nrtZeBA3ntXTw+uVC4t+tG7tiLg28AswBvhNxLN2lH0F5hPCxX4SxdJfrgfnAo8TgT+ndUELwl80s\nknn9SgjOb/tdoEm5fnsQ/HKrp6ms1aHrlwuJPx9WPD+J4B/gN4HvEpQSkixN8q7rL4GBBGWEj4D7\nog2n0/YGngJuADa0eC0J129voILg/DaSrOvXSHAeA4Ay4LQWr+/y+uVC4q8j6JDZ7lCCVn+SfBQ+\nfgw8TVDeSpp6gvoqwMHA6ghj6Q6rafoP9RvifQ0LCJL+b4Hp4b4kXb/t5/c7ms4vSddvu8+A/wOc\nQAevXy4k/jnAYJpu7LoQaGPx3FgqBL4SbvcFRrJjp2FSPAtsn0p1LE3/4ZLi4GbbFxDfa9iDoNSx\nCHig2f6kXL+2zi8p1+9AmspUfYCzCEYpxfL6JfnGroEE9bh5BMPLknB+U4APgS0E/TNXEIxaeoUY\nDSdrR8vzuxKYTDAkdz7Bf6q41sBPJigVzGPHoY1JuX6tnd83Sc71+xrwFsH5vQ38KNyflOsnSZIk\nSZIkSZIkSZIkSZIkSZIkSfH2/wFKukqXdk5eygAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10c314390>"
       ]
      }
     ],
     "prompt_number": 56
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example from the textbook: fifth order polynomial interpolation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline\n",
      "\n",
      "def evalPoly(a,xData,x):\n",
      "    n = len(xData)-1  # order of polynomial\n",
      "    p = a[n]\n",
      "    for i in range(1,n+1):\n",
      "        #print \"a[\", (n-i),\"]=\", a[n-i]\n",
      "        p = a[n-i] + (x-xData[n-i])*p\n",
      "    return p\n",
      "\n",
      "def coefficients(xData,yData):\n",
      "    n = len(xData)\n",
      "    a = np.zeros((n, n));\n",
      "    for i in range(0, n):\n",
      "        a[i, 0] = yData[i];\n",
      "    for j in range(1, n):\n",
      "        for k in range(0, n-j):\n",
      "            a[k, j] = (a[k+1, j-1] - a[k, j-1])/(xData[k+j]-xData[k])\n",
      "    return a[0, :] # return the zeroth row\n",
      "\n",
      "xAll = [0, 10, 15, 20, 22.5, 30]\n",
      "yAll = [0, 227.04, 362.78, 517.35, 602.97, 901.67]\n",
      "plt.plot(xAll,yAll,'bo',markersize=10)\n",
      "\n",
      "xData = xAll\n",
      "print \"xData=\",xData\n",
      "yData = yAll\n",
      "print \"yData=\",yData\n",
      "a = coefficients(xData, yData)\n",
      "print \"a=\",a\n",
      "\n",
      "y = evalPoly(a, xData, 16)\n",
      "print \"f(16)=\", y\n",
      "\n",
      "xs=np.linspace(0, 30, 50)\n",
      "ys=np.zeros(50)\n",
      "i=0\n",
      "for x in xs:\n",
      "    ys[i] = evalPoly(a,xData,x)\n",
      "    i=i+1\n",
      "plt.plot(xs,ys,'r--')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "xData= [0, 10, 15, 20, 22.5, 30]\n",
        "yData= [0, 227.04, 362.78, 517.35, 602.97, 901.67]\n",
        "a= [  0.00000000e+00   2.27040000e+01   2.96266667e-01   4.01666667e-03\n",
        "   6.30222222e-05   1.43407407e-06]\n",
        "f(16)= 392.070578916\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 58,
       "text": [
        "[<matplotlib.lines.Line2D at 0x10cb96110>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGhFJREFUeJzt3Xt4FfWd+PE3YECCumptUaJpUMQK1jtKtGCk6qJtUVc2\namtlq9bdxVbcPluRx/6UVlsvrUvrWqgXrNhWLMaqKGpRSp5gGzAI3gjIZY1ILEi9IHAQApzfHzOa\nEEIkJ8mZM3Per+fJc+bMmXPOZxz55JvP9zMzIEmSJEmSJEmSJEmSJEmSpBx2P7AGeK3Juv2B54Cl\nwExg3yavjQOWAUuAs5qsPyH8jGXArzoxXklSOw0BjmPHxH87cG24PBa4NVweALwMFAAlwHKgS/ja\ni8BJ4fLTwPBOi1iS1G4l7Jj4lwC9w+UDw+cQjPbHNtnuWWAwcBCwuMn6i4DfdEagkqTP1jWD9/Qm\nKP8QPn7yS6APsKrJdquAohbW14frJUkRyCTxN5UOfyRJMbFHBu9ZQ1DiWU1Qxnk3XF8PHNJku4MJ\nRvr14XLT9fUtffBhhx2WXrFiRQYhSVJeWwH0292NMxnxTwdGhcujgMebrL8I6A70BQ4nmNRdDXwE\nnEww2fvtJu/ZwYoVK0in04n9ufHGGyOPwX1z/9y/5P0Ah7UliX/WiH8qcBpwAPA2cANBF8804HKg\nDigPt60N19cCW4HRNJaBRgMPAD0JunqebUuQkqSO81mJ/+JdrD9jF+t/Fv409xLw5d0NSpLUedo7\nuas2KCsrizqETpPkfQP3L+6Svn9t1eWzN8mqdFivkiTtpi5dukAb8rkjfknKMyZ+ScozJn5JyjMm\nfknKMyZ+ScozJn5JyjMmfknKMyZ+ScozJn5JyjMmfknKMyZ+ScozmdyIRZIUsVQqRVXVXCoqqtv8\nXi/SJkkxU15+AzU1PaivL6WhYTDQC9qQzx3xS1KMpFIpamp6UFd3fcafYY1fkmKkqmou9fWl7foM\nE78kxUhFRXVY3smciV+SYmTlyk1AIV3ZlvFnmPglKUaKi3tyCffxLMOBzJphnNyVpLhIp7l225v0\n4E7O4S9k2pjpiF+S4mDbNrjqKvq99CLlB/8btQzM+KNM/JKU67ZsgX/5F1i2jK4vvMAXS3tQUnIz\nBQWzgFSbP84TuCQp16XT8OCDcPHF0L07sOOZu5Mn/wjakM9N/JIUc126dIE25HNLPZKUZ0z8kpRr\nUm2v27eFiV+SckU6DbfcAuee26lfYx+/JOWChgYYPRrmz4cZMzr1q0z8khS1deugvBy6dYOqKth7\n7079Oks9khSl9evhlFOgf3+YPr3Tkz7YzilJ0ZszB4YMyfjtbW3nNPFLUszZxy9JapWJX5KyZfNm\nePPNqKMw8UtSVqxdC2ecAXfcEXUk7Ur844BFwGvAQ0APYH/gOWApMBPYt9n2y4AlwFnt+F5JipfX\nX4eTT4bTToM774w6mownd0uAvwBHApuBPwJPAwOBfwC3A2OB/YDrgAEEvxwGAUXA80B/YHuzz3Vy\nV1KyzJgB3/kOTJgA3/pWp3xFtiZ3PwIagEKCk8AKgXeAEcCUcJspwHnh8rnA1PA9dcBy4KQMv1uS\n4uGNN+DKK+GJJzot6Wci0zN33wfuAFYCm4A/E5R4egNrwm3WhM8B+gBzm7x/FcHIX5KS64gjYPFi\n2GefqCPZQaaJ/zDgGoKSzzrgEeCSZtukaf1OwC2+Nn78+E+Xy8rKKCsryzBEScoBnZD0Kysrqays\nzPj9mdb4LwTOBK4In38bGAwMA04HVgMHAbOBLxHU+QFuDR+fBW4E5jX7XGv8ktRG2arxLyFI9D3D\nLzsDqAWeBEaF24wCHg+XpwMXAd2BvsDhwIsZfrck5Z7HHoMX45HWMi31vAI8CMwn6MxZANwD7A1M\nAy4nmMQtD7evDdfXAluB0bReBpKkeNi+HW66CSZPDpJ/DHitHknK1MaNMGoUvPMO/OlPcOCBkYTh\ntXokKRveegtOPTW4jPLs2ZEl/Uw44pekTDz5JKxYAWPGQJdoU6mXZZakPGOpR5LUKhO/JH2WzZuj\njqBDmfglqTUvvBBceuGdd6KOpMOY+CWpJek0TJwIF1wAv/kN9OkTdUQdJtMTuCQpuTZtgtGjoaYG\n/vpX6Ncv6og6lIlfkppKp2HECDjgAJg3D3r1ijqiDmc7pyQ19+abUFISeX/+7rKPX5LyjH38kqRW\nmfgl5a/aWnjggaijyDoTv6T89PDDUFYWmzp+R7KrR1J+2bIFfvhDmDEDZs6EY4+NOqKsM/FLyh/1\n9fCv/wqf/zzMnw/77ht1RJGw1CMpf2zbBiNHBnfKytOkD7ZzSlLs2c4pSWqViV9SMi1aFNwIXTsx\n8UtKlnQaJkyA00+H5cujjiYn2dUjKTk+/BAuuwzefju4wFrfvlFHlJMc8UtKhgUL4MQTg+vmv/CC\nSb8VdvVIir90Gs4/Hy6+GC68MOposs6rc0rKT+l0Xl5+AWznlJSv8jTpZ8LELyle0mnYuDHqKGLN\nxC8pPt5/P6jljx8fdSSxZuKXFA9z5sBxxwXdOj/9adTRxJp9/JJy27ZtQaKfOBEmT4avfS3qiGLP\nxC8pt91zD1RWBn36ffpEHU0i5No0uO2cknbU0ABdu0K3blFHkrPs45ekPGMfv6T42ro16gjygjV+\nSVmRSqWoqppLRUU1K1duori4JyNHljJ06GAKe/aEe++FSZOCWyJa1ulU7Sn17AvcBwwE0sB3gGXA\nH4EvAnVAOfBhuP044DJgG3A1MLOFz7TUIyVQefkN1NT0oL6+lIaGwUAhkKKgYC4DD5zFb7Y/yskH\n7AlTp8KRR0Ydbuxks9TzK+Bp4EjgaGAJcB3wHNAfmBU+BxgAXBg+DgcmtvO7JcVEKpWipqYHdXXX\n09AwjCDpAxRyakMXpr/9IK+t/xyp2bNN+lmSafL9J2AIcH/4fCuwDhgBTAnXTQHOC5fPBaYCDQR/\nCSwHTsrwuyXFSFXVXOrrS3dafxjL+T2X8F3uZfSmn1A1b2EE0eWnTBN/X2At8FtgAXAv0AvoDawJ\nt1kTPgfoA6xq8v5VQFGG3y0pRioqqsPyzo5W0I/+LOXPDKehoZSKiuoIostPmSb+PYDjCUo2xwMb\naSzrfCId/uyKxXwpD6xcuYnG8s6OUvQKlwrD7ZQNmXb1rAp/asLnFQSTt6uBA8PHg4B3w9frgUOa\nvP/gcN1Oxje5+FJZWRllZWUZhigpFxQX96Q7H7CF/VrZKkVxcc+sxRR3lZWVVFZWZvz+9nT1VAFX\nAEuB8TT+Sn8PuI3gL4B9w8cBwEMEdf0i4HmgHzuP+u3qkRJm/k238PkbfsmpvEQ9B7e4TUHBLKZP\n78Lw4cOyHF0ytLWrpz19/N8H/gB0B1YQtHN2A6YBl9PYzglQG66vJZgIHo2lHinZUikYO5bjn3iC\nb/b+OvVrWk76AEVF1Qwd+oMsBpffvGSDpI43fz5ccgmccALcdRfl/z6BmpruYR9/KY19/NUUFVUz\naNAWpk37SdRRx5bX6pEUrU2bgoR/44073Pi81TN3C1ue/NXuMfFLit7WrbCHV4TJFhO/JOUZr84p\nKXtWrw7ukKVYMfFLart0Gv7wBzj6aJg3L+po1EYW4SS1zdq18J//CYsXwzPPBBO5ihVH/JJ23/Tp\ncMwx0LcvvPSSST+mHPFL2j3pNDz6KEybBl/5StTRqB3s6pGkmMvmJRskxZgnVOUvR/xSHmrtVohF\nRdVcVryY//fHX8CBB0YdqnaDffySWtXarRB7NAzi2rpVXPbXp/h42bIow1QnMvFLeWZXt0I8nb/w\nCsfQnS0c02UKlRsbIohO2WDil/LMzrdCTHMXVzGFUXyPu7iCyby39UxvhZhgJn4pz+x8K8QuPM8Z\nHMXrPMM54TpvhZhkdvVIeSa4xWGKpsn/cc5vtpW3QkwyR/xSnhk5spSCgrmtblNQUM3IkTvPAygZ\nTPxSvli7Fr75TcrSmykqar1+H9wKcXCr2yi+TPxS0qXT8NBD8OUvQ1ERe552GoMGbaak5GYKCmYR\nlH0g6OOfRUnJzQwatMWTuBLME7ikJHv77eBKmitXwuTJMGjQpy955m5yeAcuSYHt2+HEE+H882Hs\nWOjePeqI1ElM/JIabd4MPXpEHYU6mYlfkvKM1+qR8tH8+cHoXtoNJn4pztavh6uvhm98A5YujToa\nxYSJX4qrp56CgQNhwwZYtCho15R2g5dskOJm0ya49FJ4+WX47W/hq1+NOiLFjIlfips994SvfQ0e\nfBB6ej0dtZ1dPZIUc3b1SEniQEidwMQv5aoZM+Coo+Ddd6OORAljjV/KNatWwZgx8OqrMHEifOEL\nUUekhHHEL+WKrVthwgQ49tigNfO11+DMM6OOSgnkiF/KFW+9Bc8/D3/7G/TvH3U0SjC7eiQp5uzq\nkSS1qr2JvxuwEHgyfL4/8BywFJgJ7Ntk23HAMmAJcFY7v1eKr1degR/9KOoolMfam/jHALXAJ/WZ\n6wgSf39gVvgcYABwYfg4HJjYAd8txctHH8E11wQTtsXF9ugrMu1JvgcD5wD30VhbGgFMCZenAOeF\ny+cCU4EGoA5YDpzUju+W4uOTe94eeWRwQbXaWrjySuiSa1Nsyhft6eqZAPwQ2KfJut7AmnB5Tfgc\noA8wt8l2q4Cidny3FB/33QeTJsEjj8App0QdjZTxiP/rwLsE9f1dDVvSNJaAdvW6lHyXXgo1NSZ9\n5YxMR/ynEJR1zgH2JBj1/45glH8gsBo4iOCXA0A9cEiT9x8crtvJ+PHjP10uKyujrKwswxClHOE9\nb9XBKisrqayszPj9HVFkPA34b+AbwO3Ae8BtBBO7+4aPA4CHCOr6RcDzQD92HvXbx6/4WrAANm6E\nIUOijkR5Jqo+/k+y9a3AmQTtnMPC5xB0/kwLH58BRmOpR0nxj3/Af/wHnH02rF4ddTTSZ8q1tgJH\n/MoZqVSKqqq5VFRUs3LlJoqLezJyZClDhw6msLAwuLbO3XfDj38MF10UPO63X9RhKw+1dcRv4pda\nUF5+AzU1PaivL6WhYTBQCKQoKJhLUVE1gwZtZtqWV2HdOrjzTu93q0iZ+KV2SqVSDBw4gbq663e5\nTUnJzSx6/mIKDz3UfnxFzmv1SO1UVTWX+vrSVrepry+latlbJn3Fkolfaqaiojos7wCkuYAK9mTT\nDts0NJRSUVGd/eCkDuD1+KVmVq7cBBRyHAu4k6vpxUYWcDxvcmiTrQrD7aT4ccQvNXPUAVu5l3/j\nac5hCqM4kfnNkj5AiuLinpHEJ7WXI36pqbfe4tYZk7ir6z/zpe1LWLfDlcUbFRRUM3Jk6/MAUq5y\nxC819cUvsu2ll/jf4mN2mfQBioqqGTp08C5fl3KZI36pmZ79+jFo0Gbg5rCPv5TGPv7qsI9/S3AS\nlxRDudaLZh+/sqO+HubOhQsu2OUmn3nmrpQjPIFLas2GDXD77fDrX8OYMXDDDVFHJLWbJ3BJLdm2\nDSZPhiOOgP/7P1i40KSvvGWNX/lh7FiYNw8eewxO8q6fym+WepQfNmyAXr28xIISyRq/JOUZa/zK\nXxs2wPjx8MYbUUci5TQTv+Jv61a45x7o3x+WLQtKOpJ2ycldxVc6DdOnw7hx0Lt3sHziiVFHJeU8\na/yKrzVr4Nxzg7bMs8924lZ5y8ldScozTu4qmRwQSB3GxK/c9uGHcP31UF4edSRSYpj4lZs+/hh+\n8YugU+fvf4c77og6Iikx7OpR7pk6NbjEwvHHw+zZMHBg1BFJiWLiV+7p2hUefhhOOSXqSKREsqtH\nkmLOrh7Fx8KFsH171FFIecfEr+xbtAjOPx9GjICVK6OORso7Jn5lz/Ll8O1vw7BhMGRIcF2dkpKo\no5LyjpO7ylib7kk7Z04wyh8zBiZOhL33jiZoSU7uKjPl5TdQU9OD+vpSGhoGA4VAioKCuRQVVTNo\n0GamTftJ4xsaGmD9eth//6hClhKrrZO7jvjVZqlUipqaHtTVXd/slUIaGoZRVzcMuJlUKtU48i8o\nMOlLOcIav9qsqmou9fWlO63/Jz7kx9zApUyhvr6Uqqq5EUQn6bOY+NVmFRXVYXknsDcf8SNuYjn9\nOJhVVDGUhoZSKiqqI4xS0q6Y+NVmK1duAgrZgwau5TaW048jeINSqrmc+6mjL1AYbicp11jjV5sV\nF/cEUmylJ/vxAaczm1qaX08nFW4nKddkOuI/BJgNLAJeB64O1+8PPAcsBWYC+zZ5zzhgGbAEOCvD\n71UOGDmylIKCuUAXxnFrC0kfCgqqGTly53kASdHLNPE3AP8FDAQGA1cBRwLXEST+/sCs8DnAAODC\n8HE4MLEd361s27QJXnzx06dDhw6mqKj1+n1RUTVDhw5udRtJ0cg0+a4GXg6XNwCLgSJgBDAlXD8F\nOC9cPheYSvALow5YDpyU4XcrWzZtgl/9Cvr1g0mTPl1dWFjIoEGbKSm5mYKCWUAqfCVFQcEsSkpu\nZtCgLTufxCUpJ3REjb8EOA6YB/QG1oTr14TPAfoATXv7VhH8olAuSqXg7rvh5z+Hk0+GJ58Mro3f\nxLRpP2ly5u6EZmfu/sCkL+Ww9ib+vYBHgTHA+mavpcOfXWnxtfHjx3+6XFZWRllZWbsCVAa++93g\nDlhPPw3HHrvLzQoLCxk+fBjDhw/LYnCSKisrqayszPj97blkQwHwFPAM8Mtw3RKgjKAUdBDBBPCX\naKz13xo+PgvcSPBXQlNesiEXNDQEZ9pKioVsXY+/CzAZqKUx6QNMB0aFy6OAx5usvwjoDvQFDgde\nRNFqaGh5vUlfSrRME/+pwCXA6cDC8Gc4wYj+TIJ2zmE0jvBrgWnh4zPAaFovA6kzvf8+3Hgj9O0L\n69ZFHY2kLPPqnPlk7Vr4n/+Be+6B886DceOCjh1JseatF9Wyhx+GI46ADz6Al16CyZNN+lKecsSf\nL1atgi5doMguWilp2jriN/FLUsxZ6slnCxZAeTm8/nrUkUjKYSb+uEunoaoKhg+HESOgtNQbmEtq\nlZdljrPaWrjySli9Gq67Dp54Anr0iDoqSTnOGn+crV4NlZUwciTs4e9wKV85uStJecbJ3aT58EO4\n5RZ44YWoI5GUECb+XPXOO3DttXDYYbB4MXzhC1FHJCkhTPy55t134Yor4KijgksjL1gADz4I/ftH\nHZmkhHBGMNcUFsKhh8KyZfC5z0UdjaQEcnJXkmLOyd04+PhjuPdeeOaZqCORlIdM/Nn03ntw003B\nmbVPPGEpR1IkTPzZ8NFH8L3vweGHQ10dzJoFTz0FJ50UdWSS8pCTu9lQWAh9+sCiRXDQQVFHIynP\nObkrSTHn5G5U1q+HO++E3/0u6kgkqVUm/vZ6++3gDNu+fWHOHBgwIOqIJKlVJv5MbdwIF10Exx4L\nW7dCTQ088giccELUkUlSq6zxZyqdDso6550H++wTdTSS8piXZZakPOPkbkdaujTov7/ppqgjkaQO\nY+JvLp2GmTPhnHNgyBDYbz+4/PKoo5KkDuMJXE2lUsHZtF27wjXXwJ/+BHvuGXVUktShrPE398or\ncPTR0CXX/tNIUsuc3N29bwnaMffaq/O/S5I6mZO7rdmwASZNCu5u9bOfRR2NJEUiP2r8y5bBr38d\n9N2fdhrcdReUlUUdlSRFIvmJ/6OP4KtfhW99CxYuhOLiqCOSpEjlR41/2zbo1q3jP1eSckD+1vgX\nLAg6clpi0pekT8U78W/eDL//PZSWBtfMWb486ogkKefFs9Tz4Ydw221w//1Bz/1VV8HXvw57JH/K\nQpKay/VSz3BgCbAMGNvSBldc8VOeffYvpFKpXX9K9+6wfXtw/fvnngtG+yZ9Sdot2Uz83YC7CJL/\nAOBi4MjmG02e/F+MGAEDB06gvPyGlj+psDAY8ffv34nhdrzKysqoQ+g0Sd43cP/iLun711bZTPwn\nAcuBOqABeBg4d+fNCmloOJ1D6oZwydOPs/mBB7IYYudK8v98Sd43cP/iLun711bZTPxFwNtNnq8K\n1+3g+9zJIgZyN/9O5aYhzOm1X9YClKR8kM3Ev1sN+qfwN0YzkQHUMmH7z3n4z693dlySlFey2dUz\nGBhPUOMHGAdsB25r3OSwNKzIYkiSlAgrgH5RB9GSPQiCKwG6Ay/TwuSuJClZzgbeIJjkHRdxLJIk\nSZKy5TNP7Iq5OuBVYCHwYrShdIj7gTXAa03W7Q88BywFZgL7RhBXR2lp/8YTdKItDH+G7/y2WDgE\nmA0sAl4Hrg7XJ+X47Wr/xpOM47cnMI+gVF4L3BKuj93x60ZQ+ikBCkhm7f9NggOTFEOA49gxMd4O\nXBsujwVuzXZQHail/bsR+EE04XSoA4Fjw+W9CEqvR5Kc47er/UvK8QMoDB/3AOYCX6GNxy8XLtK2\nmyd2xV6uXRepPeYAHzRbNwKYEi5PAc7LakQdq6X9g2Qcw9UEgyuADcBigvNpknL8drV/kIzjB/DJ\n9Wy6EwycP6CNxy8XEv9undgVc2ngeWA+8N2IY+ksvQnKI4SPvSOMpbN8H3gFmEwM/pTeDSUEf9nM\nI5nHr4Rg/+aGz5Ny/LoS/HJbQ2NZq03HLxcSfxburh65Uwn+BzwbuIqglJBkaZJ3XCcBfQnKCH8H\n7og2nHbbC3gUGAOsb/ZaEo7fXkAFwf5tIFnHbzvBfhwMDAVOb/b6Zx6/XEj89QQTMp84hGDUnyR/\nDx/XAo8RlLeSZg1BfRXgIODdCGPpDO/S+A/qPuJ9DAsIkv7vgMfDdUk6fp/s3+9p3L8kHb9PrANm\nACfQxuOXC4l/PnA4jSd2XQhMjzKgDlYI7B0u9wLOYsdJw6SYDowKl0fR+A8uKQ5qsnw+8T2GXQhK\nHbXAL5usT8rx29X+JeX4HUBjmaoncCZBl1Isj1+ST+zqS1CPe5mgvSwJ+zcVeAfYQjA/8x2CrqXn\niVE7WSua799lwIMELbmvEPyjimsN/CsEpYKX2bG1MSnHr6X9O5vkHL8vAwsI9u9V4Ifh+qQcP0mS\nJEmSJEmSJEmSJEmSJEmSJEmKt/8PExtp3sWrOW0AAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10ddd6dd0>"
       ]
      }
     ],
     "prompt_number": 58
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Divided difference is invariant under all permutations the data points"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# add the data points (0,8) (8,0) (3,2) (4,5) (1,12)\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline\n",
      "\n",
      "def evalPoly(a,xData,x):\n",
      "    n = len(xData)-1  # order of polynomial\n",
      "    p = a[n]\n",
      "    for i in range(1,n+1):\n",
      "        #print \"a[\", (n-i),\"]=\", a[n-i]\n",
      "        p = a[n-i] + (x-xData[n-i])*p\n",
      "    return p\n",
      "\n",
      "def coefficients(xData,yData):\n",
      "    n = len(xData)\n",
      "    a = np.zeros((n, n));\n",
      "    for i in range(0, n):\n",
      "        a[i, 0] = yData[i];\n",
      "    for j in range(1, n):\n",
      "        for k in range(0, n-j):\n",
      "            a[k, j] = (a[k+1, j-1] - a[k, j-1])/(xData[k+j]-xData[k])\n",
      "    return a[0, :] # return the zeroth row\n",
      "\n",
      "xAll = [8, 4, 0, 1, 3]\n",
      "yAll = [0, 5, 8, 12, 2]\n",
      "plt.plot(xAll,yAll,'bo',markersize=10)\n",
      "\n",
      "xData = xAll\n",
      "print \"xData=\",xData\n",
      "yData = yAll\n",
      "print \"yData=\",yData\n",
      "a = coefficients(xData, yData)\n",
      "print \"a=\",a\n",
      "\n",
      "\n",
      "xs=np.linspace(-0.2, 8.2, 50)\n",
      "ys=np.zeros(50)\n",
      "i=0\n",
      "for x in xs:\n",
      "    ys[i] = evalPoly(a,xData,x)\n",
      "    i=i+1\n",
      "plt.plot(xs,ys,'r--')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "xData= [8, 4, 0, 1, 3]\n",
        "yData= [0, 5, 8, 12, 2]\n",
        "a= [ 0.         -1.25       -0.0625      0.2172619  -0.23988095]\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "[<matplotlib.lines.Line2D at 0x1072a2990>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEACAYAAACj0I2EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHIJJREFUeJzt3Xl8lNW9x/FPgASMIEipLJE0LmUTWUQQUHGkqFFblRa9\nWq9L9WpvK1htqyxeIBFqxdai92ptFbBwtRQuVS7gBQTqiC2LKBARZLEQgVgWQRadCiHk/nEm+4RM\nMvPMeZbv+/WaV2aeeWbOz5H8cuY85/wOiIiIiIiIiIiIiIiIiIiIiIiIiIh4QjNgNbAe2AT8Mnq8\nNbAE2Aq8CbSyEp2IiNRLZvRnE2AVcBnwFPBo9PhI4EkLcYmISANlAmuAC4DNQNvo8XbRxyIi4nKN\nMEMuRzE9c4DPKz2fVu2xiIi4XEvMkMuV1EzgB1MfjohIsDRJ4nsdBt4A+gB7MUMte4D2wL7qJ/fs\n2bO0oKAgic2LiARCAdAr1hONEnzjNlTMYDkNuApYB8wD7ooevwuYWyOiggJKS0ut38aPH289Brfc\n9Fnos9Bn4f7PAuhZW0JOtIfeHpiO+cPQCPhvYFk0qc8G7gUKgVsSbEdEROqQaELfAFwU4/hBYEiC\n7y0iIvWQ6JCL54VCIdshuIY+iwr6LCo0+LM4eTL28ZUroVcvePRR2L+/wXHZ4PZ/F2kW2y6NjgeJ\niJ/s2wfPPgtvvAFr10Kjav3GI0dg82aYMQNmzoT77oOf/xzatLETr8ekpaVBLbk78D10EUmSTz6B\nESOgSxc4eBBee61mMgc44wzo1w+eew7WrYNDh6BzZ1izJvUx+4x66CKSuBkz4OGH4d/+DR56CNq3\nr9/rCwuhQwfIyHAkPD85VQ9dCV1EErd+PbRoAeedZzsS31NCFxHxCY2hi4j37NgBR4/ajsJTlNBF\npH4+/TQ17fzud/Dd78Lx46lpzweU0EUkfgUFZg75zp3Ot/WLX8Dpp8M999Q+p12qUEIXkfjs3g3f\n/jY8/zxkZzvfXpMmZp76jh3w2GPOt+cDuigqInU7ehQuvxxuuw1Gjkxt2wcOQLdusGwZdO+e2rZd\nSLNcRKThTpyAG2+ErCz4/e8hzULamDIFzjwTvve91LftMqdK6Mmshy4ifrRnj1n08/zzdpI5mAVL\nUif10EVEPETz0EVEAkAJXUTEJxJN6B2Bt4CNwIfAg9HjecBuzM5F64DcBNsRkVRx+1Boaakpvys1\nJDqG3i56Ww80B94HbsJsOXcU+M0pXqsxdBE3mjDBFNp66CHbkcRWWAgXXwwffQRf/7rtaFLOyTH0\nPZhkDvAF8BGQVdZugu8tIqm2apWpU37zzbYjqV1ODnz/+/D447YjcZ1kjqHnAL2BVdHHI4ACYCrQ\nKontiIgTjhyB2283NVSysuo+36Zx48wq0i1bbEfiKsnqRTcHwsBEYC5wFlC2WeAEoD1wb7XXlI4f\nP778QSgUcv1+fSK+dvfdkJ4OL71kO5L4PPEEbN0Kf/iD7UgcFQ6HCYfD5Y/z8/PBwZWi6cACYCHw\nTIznc4D5wIXVjmsMXcQt5s41S/rXrjUFsbxg3z7o1MkUCjvjDNvRpIyTY+hpmCGVTVRN5pX3nxoK\nbEiwHRFx0pAhsGCBd5I5wFlnwbx50LSp7UhcI9Ee+mXAcuADoKy7PQa4DegVPbYD+CGwt9pr1UMX\nEaknFecSEfEJLf0XEQkAJXSRoCopsR2BJJkSukgQvfMOXHON7SiS5+hRs4I04JTQRYLm2DG4/374\n8Y9tR5I8r78Ow4fbjsI6XRQVCZr8fDPffO5cextWJNuXX0LHjrBhg/tXuSZIs1xExNi8GS67DNat\nMwnQT+6/H845B0aPth2Jo5TQRcSUnQ2FzL6cDz5Y5+mes3o1/Ou/mnIAfvnmEYOmLYqIMXo0PPCA\n7Sic0a8fZGSYC74BpR66iPjHG2+YMfRevWxH4hgNuYiI+ISGXEREAkAJXUTEJ5TQRfxsyRJ48UXb\nUUiKKKGL+NXx4zBiBHToYDsSSREldBG/eu45s9Dm+uttR5J6L70EEybYjiLllNBF/GjvXrPn5uTJ\nvl5kU6vOnU1pg4BJNKF3BN4CNgIfAmXLz1oDS4CtwJtAqwTbEZH6GDPGbPrcpYvtSOwYMAC2b4c9\ne2xHklKJJvRi4GHgAqA/8ADQFRiFSeidgGXRxyKSCidOmIqKY8fajsSe9HSzT+qiRbYjSalEE/oe\nYH30/hfAR0AWcAMwPXp8OnBTgu2ISLyaNIFXXoGWLW1HYtd118HChbajSKlkjqHnAL2B1UBbKjaF\n3ht9LCKSOrm5pkxwgFakN0nS+zQH/gz8BDha7bnS6K2GvLy88vuhUIhQKJSkcEQk8Nq3hy1bPH9R\nOBwOEw6H4zo3Gf+l6cACYCHwTPTYZiCEGZJpj7lwWv3qjGq5iIjUk5O1XNKAqcAmKpI5wDzgruj9\nu4DgzR8SSaXDh7WnpiTcQ78MWA58QMWwymjgXWA2kA0UArcAh6q9Vj10kWQZORIOHjQLasTXVD5X\nxM927ICLLzb7aWqZv++pfK6In40ebbaUUzKP7cgRePNN21GkhHroIl62ahUMG2Zmc5x+uu1o3Gnv\nXlMKYP9+s+DI49RDF/GrMWNMESol89q1bQvnnw8rVtiOxHHqoYt4WWEhdOwIjRvbjsTdxo0z5RAm\nTbIdScLUQxfxq5wcJfN4BKQMgHroIuJ/JSXQpg1s3myGYDxMPXQRCbbGjeFXv4KTJ21H4ij10EVE\nPEQ9dBG/2L0b7rgjUBUEJX5K6CJeMnasmdXi8QqC4oxklc8VEacVFJiZGlu22I5EXEo9dBGvePRR\n+I//0E5EUisldBEvWLzYFOH64Q9tR+Jt778Pjz9uOwrHKKGLeMGuXfD0076oRWJVixYwdartKByj\naYsiEhylpfD1r5vrEVlZtqNpEE1bFBEBMzto4EBYudJ2JI5IRkKfBuwFNlQ6lgfsBtZFb7lJaEdE\nJHEDBiihn8LL1EzYpcBvgN7R26IktCMikriBA31bSjcZCf0d4PMYx7XyQSQRBQVaEeqESy6BadNs\nR+EIJ8fQRwAFwFSglYPtiPjPBx/AVVfB0aO2I/GfZs2ga1fbUTjCqZWiLwBlkz0nAE8D91Y/KS8v\nr/x+KBQiFAo5FI6Ix5QtIjrjDNuRiGXhcJhwOBzXuckaFskB5gMX1uM5TVsUiWXxYhgxAj78EDIy\nbEcjLmNj2mL7SveHUnUGjIjUpqQEHnnEbJWmZC71lIwhl5nAFUAbYBcwHggBvTCzXXYAWq8sEo9X\nX4VWreCmm2xH4n+lpWbDCx9t4aeVoiJucuwYfPaZZ1cxespdd0FuLtx2m+1I6kUrRUW8omlTJfNU\nufBC3y0wUkIXkWDy4YpRDbmISDB99RV87Wuwfz9kZtqOJm4achERqa5ZMzPs8t57tiNJGiV0EZtK\nSuA734GiItuRBNMVV8Df/247iqTRnqKWRCIRli9fxZw5K9m5859kZ5/GsGEDGDSoP5ke+vonCZoy\nBY4cgQ4dbEcSTJMm2Y4gqTSGbsEtt4xjzZqmFBUNoLi4P5AJREhPX0VW1kr69j3G7Nn+3SZLog4d\ngi5dzMbPvXvbjkY8QmPoqbJhA/zsZ3DeeWarq6ws8ws7e3b5KZFIhDVrmlJY+BjFxYMxyRwgk+Li\nwRQWPsaaNRlEIhEr/wmSQhMmmOEWJXNJEg25JNPHH5tE/tprcO655qv0kSNmy6uo5ctXUVQ0AIA+\nvMcmuvFPqg6xFBUNYPnyVeTmDk5p+JJCmzbBjBmmXotIkqiHnkxDh0JeHvTsWdFD79oV2rQpP2XO\nnJXRYRZ4hF9RRBZTuJez2VV+TnHxAObM8df8WKnm6FF49llo29Z2JOIjSugN8dZbcPBgg166c+c/\nKRtmuZVZdGMTn9KBdfRmFL8kg2NAZvQ88a1LLoHvf992FALmj+umTbajSAol9PooKYH8fLj9digs\nbNBbZGefBlSMj++hPeOYQD/epT+rmMzDQCR6nog47oMP4Ac/sB1FUmiWS7z27TOJ/MQJmDkT2rVr\n0NssWvQXbriB6AXRmpryFSfT/8a8eWkaQxdJhcOHzfDo4cOeqLyoWS6J2rkT+vSB/v1hyZIGJ3OA\nQYP6k5VV+/j4MZqRlbWSQYP6N7gNEamHli3NxAUfLDBSQo/HnDnw4INmmlmTxCYGZWZm0rfvMXJy\nJpKevoyK4ZcI6enLyMmZSN++x6suLtq/H4qLE2pXLCsthd/+1tQPEffp0cMMvXhcMoZcpgHXA/uo\n2GauNTAL+AZQCNwCHKr2Om8NuZSWQlryRqjqtVJ05EhYsQJmzdKKQq965RV45hlYvdoTX+sDZ+xY\n8/v9uPsX9J1qyCUZGepy4AtgBhUJ/Sngs+jPkcCZwKhqr/NWQrfp5El44gnTw3v1VbjyStsRSX0c\nOgTdusHcudCvn+1oJJYlS8yagIcfth1JnZxO6FBzI+jNmG3p9gLtgDDQpdprlNDra+lSuOMOGDUK\nfvIT29FIvIYPh+PH4cUXbUciPnCqhO7UStG2mGRO9KdWTyTDkCHmK3turunxXXWV7YikLn/9K7z+\nulaESkqk4qJoafTmDSUlcN99ZmaLG2Vnw7vvmuQu7jd/PvzXf8GZZ9qORALAqR562VDLHqA95oJp\nDXl5eeX3Q6EQoVDIoXDqIT/fTF9y88XH5s1tRyDx8ll5Vkm9cDhMOByO61ynxtCfAg4AkzAXQ1vh\nhYui27ebi1YbN6rGhoi4ktMLi2YCK4DOwC7gB8CTwFXAVmBw9LH7jRtn5pt7MZlv3w5bt9qOQsS7\nDh0yq8A9TEv/yxQUwDXXwLZtplKi1/zpT6YW++LF0L277WhEvOfgQcjJMYm9kXvXXNqY5eI9Bw/C\n0097M5kD3HqrWfw0ZAj83//BRRfZjiiYtm0zq3q7dbMdidRX69amDEBhodnPwIMCmdBrXaU5dCie\n3s3zttvMTua5ufC//wsDBtiOKFhKSuDuu83/ByV0byorAeDRhB64IZdA7Oe5aBHceSe8/z507Gg7\nmuD49a/NH9K333b1V3Y5hTFjTKdo3DjbkdRKQy5RlffzrKpsP8/BwEQikUjNeipekptr5qormafO\n2rVmiuK77yqZe1mPHqYYn0cF6l9e5f08a1O2n6fn5eTYjiA4vvzSDLP853/COefYjkYSMXCgmRzh\nUYFK6JX38wToxkaaUrWcqfbzlHpbuxYGDzZJXbwtO9usFPeoQCX0yvt5ZnCMN7iePrxf7Swf7+ep\nWtzOuPxyeOEF21GIBCuhV97P8z5e4kO6s4JLq53l0/08T56ESy+F//kf25GIiEMCldCHDRtAevoq\noJQf8QJP8WiNc9LTVzJsmA+n+zVqBNOmmZWwHl8NJyKxBWqWi9nPczLtCk8jnWLe4fIa55j9PH9q\nIboU6NnTFPK/+mqz2fUdd9iOSESSKFA99LL9PB9qPpyXG4WAsrHyU+zn6Tfdu5uNMkaPhilTbEfj\nTeEw/PGPtqMQpxQXwwMPmJXXHhO4hUUAXy1cyIrDEf64dHPd+3n61bZtsGwZ/Pu/247EWwoLoX9/\nsxXgt75lOxpxSvv2rl3LkYot6BrCXcW5ROryxRfmwvI992gLQL+75hpzven6621HUoPT5XNF/O/k\nSVOnpU8f84su/tajh6nA6jFK6CLxmDQJiorMfPM0m19sJSW6dIEtW2xHUW9K6FJh2zYYNcpUDZSq\nbrgBXnsNmja1HYmkQqdOntwwxumuRiFwBCgBioF+lZ5L7Rj6sWOmcL0XdyNKlcOH4aaboE0beOWV\nwCWvWssqB+liuRhHjsDKla6s62LzougOoA9wMMZzqU3os2aZhTWLF6euTS/66iszP/3AAZg7F844\nw3ZEKRGIssriC7YvirpjwHHKFDM7QU6tWTOznV2XLnDFFbBnj+2IHFe5rHJx8WAo3+akrKzyY6xZ\nk0EkErEZpkidnE7opcBS4D3AXgmzHTtg/XoznCB1a9wYnn8evvtds7LU56qXVf4an3En06uc45uy\nyuJrTif0S4HewLXAAxBjrX0qTJsGt98euDHhhKSlwdixgSgPULmscksO8SZX05mqMxxUVlm8wOla\nLv+I/twPvI65KPpO2ZN5eXnlJ4ZCIUKhUPIjKCmBl1+GhQuT/97iC2VllZtzlIVcy3IG8Ri/qHaW\nj8sqi6uFw2HC4XBc5zqZ0DOBxsBR4HTgaiC/8gmVE7pjvvjC1GW48ELn2xJPys4+jdP4jPkM4wN6\n8DCTqXnpx6dllaV2n3xivqXOmGE1jOqd3fz8/FrPdXLIpS2mN74eWA0sAN50sL3YWrY0hagkObZt\nM3OyDxywHUnSDBs2gKlpt7KTbH7EC8S6ju/bsspSu9at4c9/NquEPcLJhL4D6BW9dQd+6WBbkirn\nngtdu0LfvrBhg+1okmLQoP48m9WHu/kDpbX8Spiyyv1jPic+1aKF6RAWFdmOJG5aKSr107ixWQY/\ncSJceaWZDePxImuZmZlkD2jKN3KeID19GWW7WgWqrLLE5rEVo6q2KA23ZQvceSecfTbMmeP5Gida\nKSo13H8/9O4NP/qR7UjKqXyuOOfECVOVrk8f25HE58QJWLBAaxIkPr/+tRlymTzZdiTllNBFANat\ng/vuM7Vq5s2DjAzbEYnb7d0Lx4+7aqML20v/JYhOnjQF0dwgEoGRI02hpeHDzZoEJXOJR9u2rkrm\ndQnUJtGSQn/5ixl/nDgRbr0VGtXed3B07HrjRrjxxopZOaq2KT6mIRdxzttvwyOPmNW6Tz0Vcw9O\nx6scHjkCf/sbXHttw99DxEU0hi72lJaaGTBjxpg5va+/Xv4VNhKJcMEFkyksfKzWl+fkTGTjxp9q\nlolIlMbQxZ60NLj5Zti8GZ580uymHlW9ymEscVU5PHQInnlG9Xok8JTQJTUaN4YhQ6BJxWWbsiqH\nZ7OLxxnLUF7jHLZjqi4bMasc/uMfMHWqGaPv0QOysmD16ip/LESSZvZs+PnPbUcRF10UFWvKqhye\n5HNKSeNu/kBv1tGCo3zM+fyVy3iYZ2pWOdy2Dd56Cy65xExD7NlTs1bEOa1bw9q1tqOIixK6WGOq\nF0b4lCzGU3Hh82t8xrls5xCtiFnlcNAgcxNJhU6dzKpoD9CQi1gzbNgA0tNrjo8foA1r6Mc2OqnK\nodh39tnw+eemFLfLKaGLNYMG9Scr69S7AKnKoVjXqBGcf74Z6nM5JXSxJjMzk759j5GTM1FVDsXd\nPDLsonnoYp2qHIqbRSIRVixcyqw3Cvhk9zHr/z5tLSzKBZ7BbEM3BZhU7XkldBFxNcdXMjfAqRK6\nU7NcGgPPAUOAImANMA/4yKH2RESSKhKJsGZN0xgrmTMpLh5MYeFgYCKRSMQ13ySdGkPvB3wMFALF\nwJ+AGx1qS0Qk6ZK2kjmFnEroWcCuSo93R4+JiHhC2UrmU4m5ktkip4Zc4hocz8vLK78fCoUIhUIO\nhSMiUj9lK5krlNKEE5wgvdKxzJormZMsHA4TDofjOtephF4EVK4K3xHTS6+ickIXEXGTspXMZUl9\nKveynEFM5+5KZ8VYyZxk1Tu7+fn5tZ7r1JDLe8A3gRwgA/gXzEVRERFPqL6SeTvn0omtVc5x20pm\npxL6CWA4sBjYBMxCM1xExEOqr2TeSic6U3VxkdtWMjtZnGth9CYi4jllK5lhIkVFA9ha3DHaQ4+Q\nnr4yOg/dXSuZtVJUROQUylYyz58Z5jev/JLhd4/jezdfGriVonVRQhcRb+nc2dTi79DBWghK6CIi\nyVBaarZVtEh7ioqIJIPlZF4XJXQREZ9QQhcR8QkldBERn1BCFxGJV2kp7NpV93mWaJaLiEi8TpyA\nzEyzYXRGhpUQNMtFRCQZmjQxc9Bd2ktXQhcRqY+cHCgstB1FTEroIiL1kZMDn3xiO4qYlNBFROpD\nPXQREZ/o2tXMdnEhzXIREfEQzXIREQkApxJ6HmYP0XXRW65D7YiISJRTOxaVAr+J3kREJAWcHHJx\nd51JERGfcTKhjwAKgKlAKwfbERFJraIic3OZRIZclgDtYhx/DHgBeDz6eALwNHBv9RPz8vLK74dC\nIUKhUALhiIikyIsvmp/5+Y43FQ6HCYfDcZ2bimGRHGA+cGG145q2KCLeNG0avP02TJ+e8qZtTFts\nX+n+UGCDQ+2IiKSeS5f/OzXLZRLQCzPbZQfwQ4faERFJPZcu/9dKURGR+jp+HFq0gC+/NCV1U0gr\nRUVEkikjA667Dg4fth1JFeqhi4h4iHroIiIBoIQuIuITSugiIj6hhC4i4hNK6CIiDXH8OCxaZDuK\nKjTLRUSkISzNRdcsFxGRZMvIgLPOclXVRSV0EZGGclkJACV0EZGGUkIXEfEJlyX01FaVERHxk4ED\nYc8e21GU0ywXEREP0SwXEZEAUEIXEfGJRBL6zcBGoAS4qNpzo4FtwGbg6gTaEBGROCVyUXQDZr/Q\n31c73g34l+jPLGAp0Ak4mUBbIiJSh0R66JuBrTGO3wjMBIqBQuBjoF8C7YiIuNfKlfDhh7ajAJyZ\nttgBWFXp8W5MT11ExH8WLIBmzaB7d9uR1JnQlwDtYhwfA8yvRzsx5yfm5eWV3w+FQoRCoXq8pYiI\nC3zjG7B6tWNvHw6HCYfDcZ2bjHnobwE/A9ZGH4+K/nwy+nMRMB6o/l+seegi4n1vvglPPQVLl6ak\nuVTMQ6/85vOAW4EM4Bzgm8C7SWpHRMRdXLT8P5GEPhTYBfQH3gAWRo9vAmZHfy4EfkwtQy4iIp6X\nnQ333287CkBL/0VEPEVL/0VEAkAJXUTEJ5TQRUR8QgldRMQnlNBFRHxCCV1ExCeU0EVEfEIJXUTE\nJ5TQRUR8QgldRMQnlNBFRHxCCV1ExCeU0EVEfEIJXUTEJ5TQRUR8IpGEfjOwESgBLqp0PAf4J7Au\nevttAm2IiEicEknoGzC7Fi2P8dzHQO/o7ccJtOG4eDdfDQJ9FhX0WVTQZ1HB7Z9FIgl9M7A1WYHY\n4vb/Qamkz6KCPosK+iwquP2zcGoM/RzMcEsYuMyhNkREpJImdTy/BGgX4/gYYH4tr/kU6Ah8jhlb\nnwtcABxtYIwiIhKHZGwS/RbwM2BtPZ9fD/RMQvsiIkFSAPSK9URdPfR4Vf7D0AbTOy8BzgW+CWyP\n8ZqYAYmISOoNBXZhpijuARZGj38P+BAzhv4+cL2V6ERERERExFtyMVMvtwEjLcdiW0fMtY6NmG9X\nD9oNx7rGmG+YtV34D4pWwBzgI2AT0N9uOFaNxvx+bAD+CDS1G45U1hiz+CkHSMdcoO1qMyDL2lFx\nTaM5sIVgfx4/BV4F5tkOxLLpwD3R+02AlhZjsSkHcx2wLInPAu6yFs0pBLWWSz9MQi8EioE/ATfa\nDMiyPZg/agBfYHpkHeyFY9XZwHXAFJIzC8yrWgKXA9Oij08Ah+2FY9URTJ7IxPxhywSKrEZUi6Am\n9CzMBd0yu6PHxPRGegOrLcdhy2TgEeCk7UAsOwfYD7yMmXL8EiaRBdFB4GlgJ2adzSFgqdWIahHU\nhF5qOwCXao4ZM/0JpqceNN8G9mHGz4PcOwfTE70IU1zvIuBLYJTViOw5D3gI09npgPk9ud1mQLUJ\nakIvwlwILNMR00sPsnTgz8ArmNW9QTQQuAHYAcwEBgMzrEZkz+7obU308RyqVlUNkouBFcABzNDT\na5h/K+ISTYC/Y/7iZqCLommYxDXZdiAucgWa5bIc6BS9nwdMsheKVT0xs79Ow/yuTAcesBqR1HAt\nZjbHx5gpSUF2GWbMeD0VdexzrUZk3xVolktPTA+9ANMrDeosF4BHqZi2OB3zjVZERERERERERERE\nRERERERERERERERERERExL7/B9P3KAGgjyYOAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x107195c10>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}