{
 "metadata": {
  "name": "trapezoid_rule_solution.v3"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Basic numerical integration: the trapezoid rule\n",
      "===============================================\n",
      "\n",
      "**Illustrates**: basic array slicing, functions as first class objects.\n",
      "\n",
      "In this exercise, you are tasked with implementing the simple trapezoid rule\n",
      "formula for numerical integration. If we want to compute the definite integral\n",
      "\n",
      "$$\n",
      "     \\int_{a}^{b}f(x)dx\n",
      "$$\n",
      "\n",
      "we can partition the integration interval $[a,b]$ into smaller subintervals,\n",
      "and approximate the area under the curve for each subinterval by the area of\n",
      "the trapezoid created by linearly interpolating between the two function values\n",
      "at each end of the subinterval:\n",
      "\n",
      "<img src=\"http://upload.wikimedia.org/wikipedia/commons/thumb/d/dd/Trapezoidal_rule_illustration.png/316px-Trapezoidal_rule_illustration.png\" /img>\n",
      "\n",
      "The blue line represents the function $f(x)$ and the red line\n",
      "is the linear interpolation.  By subdividing the interval $[a,b]$, the area under $f(x)$ can thus be approximated as the sum of the areas of all\n",
      "the resulting trapezoids. \n",
      "\n",
      "If we denote by $x_{i}$ ($i=0,\\ldots,n,$ with $x_{0}=a$ and\n",
      "$x_{n}=b$) the abscissas where the function is sampled, then\n",
      "\n",
      "$$\n",
      "   \\int_{a}^{b}f(x)dx\\approx\\frac{1}{2}\\sum_{i=1}^{n}\\left(x_{i}-x_{i-1}\\right)\\left(f(x_{i})+f(x_{i-1})\\right).\n",
      "$$\n",
      "\n",
      "The common case of using equally spaced abscissas with spacing $h=(b-a)/n$ reads simply\n",
      "\n",
      "$$\n",
      "   \\int_{a}^{b}f(x)dx\\approx\\frac{h}{2}\\sum_{i=1}^{n}\\left(f(x_{i})+f(x_{i-1})\\right).\n",
      "$$\n",
      "\n",
      "One frequently receives the function values already precomputed, $y_{i}=f(x_{i}),$\n",
      "so the equation above becomes\n",
      "\n",
      "$$\n",
      "   \\int_{a}^{b}f(x)dx\\approx\\frac{1}{2}\\sum_{i=1}^{n}\\left(x_{i}-x_{i-1}\\right)\\left(y_{i}+y_{i-1}\\right).\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Exercises\n",
      "\n",
      "### 1\n",
      "\n",
      "Write a function `trapz(x, y)`, that applies the trapezoid formula to pre-computed values, \n",
      "where `x` and `y` are 1-d arrays."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def trapz(x, y):\n",
      "    return 0.5*np.sum((x[1:]-x[:-1])*(y[1:]+y[:-1]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### 2 \n",
      "\n",
      "Write a function  `trapzf(f, a, b, npts=100)` that accepts a function `f`, the endpoints `a`\n",
      "and `b` and the number of samples to take `npts`.  Sample the function uniformly at these\n",
      "points and return the value of the integral."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def trapzf(f, a, b, npts=100):\n",
      "    x = np.linspace(a, b, npts)\n",
      "    y = f(x)\n",
      "    return trapz(x, y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### 3\n",
      "\n",
      "Verify that both functions above are correct by showing that they produces correct values \n",
      "for a simple integral such as $\\int_0^3 x^2$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "exact = 9.0\n",
      "x = linspace(0, 3, 50)\n",
      "y = x**2\n",
      "\n",
      "print exact\n",
      "print trapz(x, y)\n",
      "\n",
      "def f(x): return x**2\n",
      "print trapzf(f, 0, 3, 50)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "9.0\n",
        "9.00187421908\n",
        "9.00187421908"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### 4\n",
      "\n",
      "Repeat the integration for several values of `npts`, and plot the error as a function of `npts` \n",
      "for the integral in #3."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "npts = [5, 10, 20, 50, 100, 200]\n",
      "err = []\n",
      "for n in npts:\n",
      "    err.append(trapzf(f, 0, 3, n)-exact)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.semilogy(npts, np.abs(err))\n",
      "plt.title(r'Trapezoid approximation to $\\int_0^3 x^2$')\n",
      "plt.xlabel('npts')\n",
      "plt.ylabel('Error')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 25,
       "text": [
        "&lt;matplotlib.text.Text at 0x472cc50&gt;"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEeCAYAAACZlyICAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVWX+B/DPBVcUBTdEURZDAUFAQxQFr1Zo4jaZpTWk\n4EYZbvOz5ufUC3KaMWo0jZlxSc2yySkbHdFflGQBYogouKQyKopUIpIYmygIz++PM964stwLdzl3\n+bxfr/uKu3DO9xyIj895nvM8CiGEABERkQY2chdARETmgYFBRERaYWAQEZFWGBhERKQVBgYREWmF\ngUFERFphYBARkVbayV0AEVFD165dw/Hjx3H+/HlMmTIFI0aMkLsk+i+2MIjIpBw9ehQ9e/bE0KFD\ncfHiRbnLoQYYGERkUp577jm4u7vjxIkTmDlzptzlUAMMDCKSXUVFBR577DGEhYXh3r17cHd3x4wZ\nMxAfHy93adQAA4OIZLd+/Xp4eXnh9u3beO6553D+/Hl07NgR165dk7s0akDByQeJSE719fVwcXHB\noUOH0L9/f1y6dAk3btxAZmYmfvvb32Lo0KFyl0j/xcAgIlllZmZi1qxZ+PHHH+UuhTTgJSkiktW3\n336LMWPGyF0GaYGBQUSySktLQ1BQkNxlkBYYGEQkm7q6OmRmZmL48OFyl0Ja4J3eRCSbU6dOobKy\nEv7+/nrZ3qVLl/D999/jzJkzmDp1KoNIz9jCINn5+voiPT29yffmzZuH119/3cgV6VdLx2cJ+9NF\nZmYmnJyc0LNnT71s7+DBg+jfvz9WrlyJv/zlL3rZJv2KgWHmunbtCnt7e9jb28PGxgZ2dnaq57t3\n75a7PK18//33CAsLa/I9hUIBhUJh5Ir0q6Xj05Wbmxu++eYbo+yvqX3pKisrC76+vnrb3ooVKzBy\n5Ej88MMPcHd319t2ScJLUmausrJS9bW7uzu2b9+OCRMmNPnZ+/fvo1078/uRyz3y25TPm0KhMNr5\nMcS+jh07hilTpuh1mwCwb98+/OEPf9D7dq0dWxgWzs3NDX//+98REhICBwcH1NXV4a233sIjjzyC\nnj174vnnn8eRI0fUPv/Xv/4VI0eOxKBBg7B582bU1taq3r99+zY2bNiAoUOH4sknn8ShQ4cAAJ9+\n+qmqZWNvb4+OHTti/PjxAICffvoJr732Gtzc3BAVFYXc3NxGNR4+fBiANFPpokWL0LdvXyxcuBD3\n799v8fjaeiyajvPh81ZfX9/sceTn56Nnz56q59evX0fv3r1Vl4UaHt+D55s2bcLo0aPh5OSEV155\nBVVVVXjmmWfg7OyM5cuXo6KiQuMxRkZGorCwEFOnToW9vb3qEkzD/bV07t3c3LB582aMHj0aAwcO\nRHx8vNo5aKipfV2/fr3Fn6smpaWlyM/Ph5+fX6u+T5OkpCTExsaisLBQr9slAIIshpubmzh8+HCj\n13x8fER6erq4e/euEEKIPXv2iKKiInHnzh2xfv164eLiovq8q6urGDJkiDhy5Ig4deqUCAwMFJs3\nb1a9/5vf/EYsXbpU3LhxQ6Snp4t+/fqJS5cuqe2zvLxceHt7i61btwohhAgLCxMvv/yyuHnzpti+\nfbvo1q2buHPnTpN1jxgxQvzud78TJSUl4p133hEdOnQQr7/+erPH3NZj0XScTZ23po6jurpaCCHE\n+++/L3x8fMSdO3dEeHi4WLVqVbM/Fzc3NxEYGChyc3PF6dOnRbdu3URQUJBISkoS169fF8HBweKj\njz7S6hib+5k/eK2lmt3c3IS/v784fvy4uHjxonBzcxNff/11s+f64X1p+rlqcujQIaFQKER2drbW\n36PJ3r17xaOPPioef/xx8cc//lFv2yUJA8OCNPfHY82aNc1+T319vRgwYIA4efKk6vMN/0Bv2bJF\nTJkyRQghBYGzs7PaH4Xly5eLt99+W/W8rq5OREREiJdeekkIIURJSYno3LmzqKysVH1mzJgxYu/e\nvY3qvnHjhujUqZPqD5oQQgwYMKDFwGjrsbT03oP3G543bY5j2rRpwtfXV/j7+4uamppGx9fw+fr1\n61XPn3jiCfHUU0+pnv/pT38Sc+fO1foYmwsMTTW7ubmJdevWqd5bvHixePXVV5vc78P70uZ8aLJ2\n7Vpha2ur9vPW1v79+8XBgwfFK6+8IrZt2yaef/55ceHChVZvh1qHl6SsQHBwsNrzpKQkPPXUU+jX\nrx969OiBoqIinD59WvV+QECA6uvAwEBkZmYCADIyMlBSUoJ+/frB0dERjo6O2L59OzIyMlSf/8Mf\n/oCqqiq89957AKRr1B4eHujSpYvqM48++qjapaMHjh8/jkceeQSdOnVSvaZpWGRbj0XTe4D6edPm\nOBYsWIBz584hNjYW7du3b7HuhsNInZyc1J736dMHP/30k9bH2Jzmam7482p4DpydndX225ZtN/Vz\nbU5ubi4GDx6s9vPWRmFhIXx8fBAREYFDhw5h5syZmD17NgYOHNiq7VDrMTCsQMMO26qqKixcuBBz\n585FXl4eSktL0b9/f7XOzIbXonNychASEgIAGD16NHr37o3i4mLcvn0bt2/fRnl5Ofbv3w8A+Oc/\n/4lPP/0Un3/+OWxtbQEAo0aNwpUrV1BVVaXaZnZ2dpOjeEaOHInLly+jurpabf/N0eVYNL338HnT\ndByVlZVYvnw5FixYgLi4ONy+fbvZupsimulM1nSMtra2zX7v6NGjm6w5NDS0VTU80HBfrfm5Nuf0\n6dNtWk1v4MCBeOSRR1BcXIzu3bvDwcEBU6ZMgZ2dXau3Ra3DwLAyFRUVqKyshLOzM+rr67F27Vpc\nv35d9b4QAv/6179w9OhRnDlzBlu3blWNYnFwcMDYsWOxevVqXLt2DXV1dfj+++9x4sQJ5ObmIjY2\nFvv27VMbU9+rVy8EBQVh9erVuHnzJnbu3Ilz585h4sSJjWpzcnLC0KFDERcXh5KSEqxfvx7FxcUG\nOZaW3muKpuNYtmwZRo4cia1btyIiIgIxMTFa/kRapukYR4wYgZMnTzb5vT179tT63Guj4b5a83Nt\nSnV1NS5fvtymwMjLy8OpU6fwxRdfqALqiy++aPV2qPUYGFamb9++WLt2LSIjI+Hv74+amhqMHTtW\n9b5CocCSJUuwcuVKzJgxA/Pnz8e8efNU72/evBmurq54+umn0bt3byxatAhlZWVISkrCL7/8grFj\nx6pGSkVERAAA/vGPf8DOzg5BQUFITU3F4cOH0blz5ybr27NnD0pLS+Hr64u8vDw8++yzej2WqKgo\nrY6zKc0dx/79+3Ho0CFs2rQJgLS2Q05OTqvug2l4r0nDe080HWNMTAwOHjyIHj16YP369VrX3FwN\nLd3z8vC+WrPth124cAH19fVtmkPq0KFD2LdvH+rr63H37l0cOHAA/fv3b/V2qPVManrze/fu4X//\n939RXV2N6dOnY9KkSXKXZHU03cthTlo6Fks6TnO0a9cuREdHo7y8XOuQIfmZVAvj6NGjCAoKwqZN\nm7B37165yyEiA7lw4QICAgIYFmbG4IERHR0NJyenRjfnpKenw9vbG56enkhMTAQAnD17FoMGDQIA\ntY5PIrIsZ8+eVbu0RubB4IERFRWFL7/8stHry5Ytw5YtW/D111/jb3/7G37++WcMGzYMV65cAQCO\neJDJ1atXLeYyTUvHYknHaY7Onj3LRZPMkMEDIzQ0FI6OjmqvlZWVAQDCwsLg6uqK8PBwZGVlISQk\nBCdOnEBsbCyeeuopQ5dGRDIoLS1FYWEhA8MMyTKjWnZ2Nry8vFTPfXx8cOzYMURERHBKYiILl52d\nDW9vbzg7O8tdCrWSaU7BqYG5T3dNRPz/WC66DIyVZZRUUFAQ8vLyVM/PnTuHUaNGtWobQpoHiw89\nPOLi4mSvwZIePJ/NP+rq6uDu7o6MjAwIIVBeXo4NGzZg69atOHHiBM+ngR+6kiUwunfvDkAaKVVQ\nUICUlJRG8x1pEh8fj9TUVANUR0T6duTIEcyZMwfr1q2Dm5ubqv/io48+wujRoxEVFdXkjYekH6mp\nqYiPj9d5OwYPjDlz5iAkJAQXL17EgAED8MEHHwAANmzYgMWLF+Pxxx/HSy+9hF69erVqu/Hx8VAq\nlQaomIj0zc/PD3l5ecjMzMSHH36oej0vLw/Ozs5o164dSktLZazQsimVSr0EhsH7MJqbHmHcuHG4\ncOGCoXdPWmDw6hfPZ2MODg5NLrDUuXNn1USVzfVp8HyaDpO607s1eElKf/g/pH7xfGpv6NChKC4u\nxt27d5u9ysDzqTt9XZIyqbmktGXMdYyJyHAqKiqQkZGBoqIiDB8+XG19DtI/Xf92MjCIiKyErn87\neUmKiMjC8ZKU+ZVNRCQrq21hEBGRcTEwiIhIK2YbGOzDICLSDvswzK9sIiJZsQ+DiIiMgoFBRERa\nMdvAYB8GEZF22IdhfmUTEcmKfRhERGQUDAwiItIKA4OIiLRitoHBTm8iIu2w09v8yiYikhU7vYmI\nyCgYGEREpBUGBhERaYWBQUREWmFgEBGRVsw2MDislohIOxxW20LZCQlASAgQGmrEooiITByH1Tah\nrAw4fFjuKoiILItFBsaoUcCxY3JXQURkWSzyklRxMeDlBdy6BdhYZCQSEbUeL0k1wckJcHAALl6U\nuxIiIsthkYEB8LIUEZG+WXRgZGXJXQURkeWw6MBgC4OISH8sstMbAO7dA3r0AG7eBLp0MVJhREQm\nzGo7vTXd6d2xIzBsGHDihPFqIiIyRbzTW4uyV6wA+vYFXn3VCEUREZk4q21haIP9GERE+mMVgWF+\nbSgiItNj0YExcKD038JCeesgIrIEFh0YCgUvSxER6YtFBwbAwCAi0hcGBhERacWih9UCQFUV0KcP\nUFoq3ZtBRGStOKxWgy5dgMGDgVOn5K6EiMi8WXxgALwsRUSkDyYXGFevXsWCBQswa9YsvW2TgUFE\npDuTCwx3d3ds27ZNr9tkYBAR6c5ggREdHQ0nJyf4+fmpvZ6eng5vb294enoiMTHRULtX4+kJlJUB\nN24YZXdERBbJYIERFRWFL7/8stHry5Ytw5YtW/D111/jb3/7G37++Wfs2rULK1aswPXr1w1Si40N\nEBzMBZWIiHRhsMAIDQ2Fo6Oj2mtlZWUAgLCwMLi6uiI8PBxZWVmIjIzEu+++i379+qG0tBQxMTE4\ndeoUEhIS9FYPL0sREemmnTF3lp2dDS8vL9VzHx8fHDt2DBEREarXevTogc2bN2vcVsO53ZVKJZRK\nZYufHzUKeOutVpdMRGS2UlNTW1w3qLWMGhj61NrFQEaOlBZTun8faGe2R01EpL2H/zH9xhtv6LQ9\no46SCgoKQl5enur5uXPnMGrUKKPs29ERcHEBzp0zyu6IiCyOUQOje/fuAKSRUgUFBUhJSUFwcHCb\ntqVpidamsB+DiKyRyS/ROmfOHKSlpeHWrVvo06cP1qxZg6ioKKSlpSEmJga1tbVYunQpli5d2upt\nt3U+lM2bpZFSH3zQ6m8lIjJ7us4lZbaTD8bFxWnV2d3QqVPAnDnAhQuGq42IyNQ86Px+4403rDMw\n2lL2/ftSX0ZhofRfIiJrwtlqW6FdO2DECCA7W+5KiIjMj1UFBsCObyKitjLbwGjLKCmAgUFE1sfk\nR0kZki7X4YqKAF9f4OefAYVCz4UREZkw9mG0krMzYG8PXLokdyVERObFbAOjrZekAF6WIiLrwktS\nOpS9YQPwn/8AmzbpsSgiIhNntTfu6VJ2Xh4wbhxw9SpgZ6fHwoiITBj7MNrAywsICQG2b5e7EiIi\n82GVLQwAOH4cePpp4PJloEMHPRVGRGTCrLaFoUunNyCtjzF4MPDJJ/qriYjIFLHTWw9lHz4MLFki\nrZFha6uHwoiITJjVtjD0YcIEoFs34N//lrsSIiLTZ9WBoVAAq1cDa9cC5tfOIiIyLqsODACYNg2o\nrgZSUuSuhIjItJltYOja6f2AjQ3w+99LrQwiIkvETm89ll1bC3h6Art3A6NH622zREQmhZ3eetC+\nPfDKK2xlEBG1hC2M/6quBjw8gEOHAD8/vW6aiMgksIWhJ507A8uXA2+9JXclRESmiS2MBsrLpVZG\nVhYwaJDeN09EJCu2MPSoWzcgJgZ45x25KyEiMj1mGxj6Glb7sGXLgE8/Ba5f1/umiYhkwWG1Bix7\n2TJpBlu2NIjIknABJQMoLAQCAqSpz3v0MNhuiIiMin0YBjBwIDBjBvDXv8pdCRGR6WALoxn/+Q8Q\nGgpcuQJ07WrQXRERGQVbGAYyZIi07vf778tdCRGRaWALowU5OdJstvn5QMeOBt8dEZFBsYVhQMOH\nA76+wK5dcldCRCQ/tjA0SEsDFiwA8vK4jCsRmTerbWEY6sa9h4WFAX36AJ9/bvBdEREZBG/cM2LZ\nBw8Cr70G5OZKy7oSEZkjq21hGFNEhLTmd3Ky3JUQEcmHgaEFhYLLuBIRMTC0NGsWUFQEHDkidyVE\nRPJoMTDq6+vx3XffGasWk9auHZdxJSLrprHTOyAgAKdOnTJWPVoxdqf3A/fuSQssHTwIBAYaffdE\nRDoxeKf31KlT8d5776G8vLzNO7EUHTsCK1dyGVcisk4aWxhdu3bFnTt3YGNjg86dO0vfpFDIGiBy\ntTAAoLIScHcHjh4FBg+WpQQiojbhehgyiI8HfvwR2LZNthKIiFrNKIFx4cIFJCUlQaFQYNq0afDy\n8mrzDvVB7sC4dQvw9ATOnAFcXGQrg4ioVQzeh7Ft2zbMmzcPNjbSR6OiorDNgP+03r9/PxYtWoTo\n6GgcP37cYPvRRc+ewEsvAc8+K4UHEZE10NjCGDNmDA4ePAhHR0cAwO3btxEREWHw4bY3b95EXFwc\nNm3a1Og9uVsYAFBfL93Md+AA8OWXgKurrOUQEWlk8BaGg4MDbjX4Z3RpaSkcHBw0bjg6OhpOTk7w\n8/NTez09PR3e3t7w9PREYmJis9+fkJCAxYsXa9yPXGxsgLffBmJigDFjABMbeUxEpHcaWxiHDx/G\n4sWL4e3tDQDIy8vDli1bMGHChBY3fOTIEXTt2hUvvPACzp49q3o9MDAQGzduhKurKyZOnIiMjAwk\nJycjJycHq1atgrOzM1599VVMnDgRjz32WNNFm0ALo6E9e4AlS4Ddu4FmSiYikp2ufzvbtfRmfX09\nbt26hYsXL+LYsWNQKBQIDg5W9We0JDQ0FAUFBWqvlZWVAQDCwsIAAOHh4cjKykJkZCQiIyMBAO+9\n9x6++eYbVFRU4PLlyybdynhg1ixpCvRnngHWrweef17uioiI9K/FwLCxscFbb72FWbNmISQkROed\nZWdnq42w8vHxwbFjxxAREaF6benSpVi6dKnO+zK2ceOAw4eByZOB69eB//kfToVORJalxcAAgBkz\nZuCVV17B3Llz0a9fP9XrPXr0MGhhmjRcDESpVEKpVMpWywO+vsB33wFPPindp7F+PVfpIyL5pKam\n6nWhOY19GG5ublA89E9lhUKBK1euaNx4QUEBpk6dqurDKCsrg1KpRG5uLgAgNjYWkyZNUmthaFW0\nifVhPOyXX4AZM4DevaX1wDt1krsiIiIDj5Kqr69HQkICrl69qvbQJiya0r17dwDSSKmCggKkpKQg\nODi4Tdsy1hKtbeHgIA21VSiA8HDg9m25KyIia2a0JVqHDx+OkydPNmplaDJnzhykpaXh1q1b6NOn\nD9asWYOoqCikpaUhJiYGtbW1be6vMPUWxgP19cDvfgccOiQFyIABcldERNbM4FODrFmzBhUVFSbV\nh2EugfHAunXAhg3AF18AD92WQkRkNAYPjKb6MADg6tWrbd6prhQKBeLi4kyms1sbu3cDy5YBn30G\nmEnJRGQhHnR+v/HGG5yt1lx88w0wezaQmCjNQ0VEZEwG6/R+++23VV/v2bNH7b3Vq1e3eYfWbMIE\nICVF6td49125qyEiap1mA2P37t2qr//85z+rvZecnGy4irRkyqOkWuLvLy2+9P77UnDU18tdERFZ\nOoOPkgoMDFTdL9Hw66aeG5u5XpJqqLQUmDZNGjm1c6e0/CsRkSEZfLZaMowePaTLUzU1wKRJwH+n\n2SIiMlnNtjBsbW1hZ2cHAKiurlat5/3g+f37941TYRMsoYXxQF2dNHoqPR1ITgb695e7IiKyVAab\nrbaurq7NGzWG+Ph4sxpW2xxbW2nU1NtvAyEh0r0aQ4fKXRURWRJ9zSnFYbUmZNcuaZbbzz8HQkPl\nroaILA37MCxIZCTw8cfAzJnAv/4ldzVEROo0Tm9OxvXEE8BXXwFTpkjrasTGyl0REZHEbAPDUvow\nmhIYCGRkSKOnfvwRWLtWWkOciKgt2IdhfmW32s8/A1OnAoMGATt2AB06yF0REZkz9mFYsF69pGVf\nKyqAiAigvFzuiojImjEwTJydndQBPmgQEBYGFBXJXRERWSsGhhlo1w7YtAl4+mnpXo28PLkrIiJr\nxE5vM6FQAK+9Jt0JrlQCe/dK4UFEpAk7vc2vbL1JTgZeeEGa8XbGDLmrISJzYbCpQch0PfmkFBrT\npkl9Gi++KHdFRGQN2MIwY/n50r0azzwDvPmmdNmKiKg5Bl/T2xQxMH5186Z0V7iPj3SJqn17uSsi\nIlPF+zCsXJ8+wLffAiUl0k1+lZVyV0REloqBYQG6dAH27wdcXKQRVMXFcldERJbIbAPDXNf0NpR2\n7aRLUlOmSMNtL16UuyIiMhUGX9PblLEPo2Xvvw+8/rrU6ggOlrsaIjIV7MOgRhYuBLZtk1obBw/K\nXQ0RWQoGhoV6EBYLFkgtDiIiXfGSlIW7eFG60S8yEoiL470aRNaM92GQRsXFwOTJ0sJMmzdLHeRE\nZH3Yh0EaOTkBqanS6n3TpwNVVXJXRETmiIFhJeztgQMHgN69gfHjpRv9iIhag4FhRdq3Bz74AAgP\nl+7VyM+XuyIiMie8mm1lFApposL+/YHQUCApCXj0UbmrIiJzYLYtDN7prZsXXwT+/ndpBNXu3cD9\n+3JXRESGwju9za9sk5SZCaxcCfzwAxAdDcyfD7i6yl0VERkCR0mRTkaPlkLjiy+A27eB4cOlIbj7\n9gG1tXJXR0SmhC0MUnPnDvD559Ld4fn5wLx50t3iHh5yV0ZEumILg/TKzk5aL/zIEeDrr4HqamkC\nw/BwYM8eoKZG7gqJSC5sYZBGd+8Ce/dKrY7z54G5c6VWx+DBcldGRK3BFgYZXKdOwHPPSSv7HTki\nvRYaKt0AuHu3FChEZPnYwqA2qamR1tvYuhU4dUqa3HDhQsDbW+7KiKg5bGGQLDp0AGbNAlJSgGPH\npFbIhAlSy2PXLqnvg4gsC1sYpDe1tdIaHFu3AsePA88/L7U6/PzkroyIAAuc3jwvLw8bN25ETU0N\nIiIi8NRTTzX6DAPD9F27BmzfDuzYAQwYIAXHs88CXbrIXRmR9bK4wHigpqYGc+fOxe7duxu9x8Aw\nH/fvA8nJUqvj6FFg9mwpPAID5a6MyPqYbB9GdHQ0nJyc4PfQ9Yj09HR4e3vD09MTiYmJTX5vUlIS\nxo8fj2eeecZQ5ZGRtGsHTJ0qTa1+5gzQt6+0JkdQkBQiFRVyV0hE2jJYC+PIkSPo2rUrXnjhBZw9\ne1b1emBgIDZu3AhXV1dMnDgRGRkZSE5ORk5ODlatWoV+/fqpPjtt2jQkJSU1LpotDLNWVwccOiQF\nRmqq1Hm+cKE0ay6XkCUyHF3/dhpsevPQ0FAUFBSovVZWVgYACAsLAwCEh4cjKysLkZGRiIyMBACk\npaVh7969EEJg1qxZhiqPZGRrK82S++STQFGRtEbHs88C3btLwfH889LXRGRajLoeRnZ2Nry8vFTP\nfXx8cOzYMURERKheGzduHMaNG2fMskhGzs7A6tXA738PHD4s3U2+ejXwm98AixYBo0ax1UFkKsx2\nAaWGc7srlUoolUrZaiHd2dgATzwhPW7eBD78UJrTqlMnqdURGQk4OspdJZF5SU1N1eu6QQYdJVVQ\nUICpU6eq+jDKysqgVCqRm5sLAIiNjcWkSZPUWhjaYB+GdRACSEuT+jq++ELqPF+0CBg7lq0OorYw\n2VFSTen+3wvT6enpKCgoQEpKCoKDg9u0La64Z/kUCkCpBD75BLh8WVqrY9EiwMcHWL8e+PlnuSsk\nMg8mv+LenDlzkJaWhlu3bqFPnz5Ys2YNoqKikJaWhpiYGNTW1mLp0qVYunRpq7fNFob1EkK6n2Pr\nVmk98ieflEJEqWSrg0gTi71xryUMDAKkFQI//lgKj7t3pb6OefOAPn3krozINJnVJSl94iUpcnQE\nYmOlGwJ37QLy8oAhQ36dFLG+Xu4KiUyDyV+SMiS2MKg5ZWVSn8fWrdLXCxYAUVHS8F0ia8dLUkRN\nEAI4eVIKjj17pD6ORYukpWZtbeWujkgevCRF1ASFQppqZOtWoLAQmDwZiIsDPDyANWuAH3+Uu0Ii\n4+ElKfMrm0zAqVPS3eS7dwNjxkgd5ZMnS5MkElk6XpIiaoOqKuCzz6TwKCwEoqOB+fMBV1e5KyMy\nHKu9JEWkiy5dpM7w774DvvxS6iAfMUK6r2PvXmn1QCJSZ7aBwT4M0hdfX2DjRuCHH6SZcjdsAAYO\nlCZBvHJF7uqIdMc+DPMrm8zIhQvAtm3ARx8BAQFSX8eMGUCHDnJXRtR27MMgMqB794B9+6TRVufO\nSTPoLlwIDB4sd2VErcc+DCID6thRWof8m2+AjAxpGvawsF8nRbx7V+4KiYzHbAODfRhkbJ6eQEKC\nNKrq5ZelNTsGDABWrADOn5e7OqLmsQ/D/MomC3TlCrB9u7TMrIeHdDf5lClAjx5yV0bUGPswiExA\nbS3wf/8ndZSnpwM9e0rrdwwfLg3XHT6cs+iS/BgYRCamvl5a8OnkSSAn59dH166NQ6RfP7mrJWvC\nwCAyA0IAV6+qh8jJk0D79uoBMny41C/CxaDIEBgYRGZKCOlmwYdDpK6ucYi4uzNESHdWGxhxcXFQ\nKpVQKpVyl0OkN0IARUWNQ6SqqnGIPPKINMyXSJPU1FSkpqbijTfesM7AMMOyidqsuFi9P+TkSWmJ\n2oCAX0O/sx16AAAKKElEQVRkxAjphkKu90HNsdoWhhmWTaRXt241DpEbNwB/f/UQ8fbm9O0kYWAQ\nkcovvwC5ub8GSE6O1E/i66seIkOHcl4sa8TAIKIWVVRIC0c1DJErV6SWR8MQ8fMDOnWSu1oyJAYG\nEbXanTvA6dPqIXLxotQH0vBeEX9/wM5O7mpJX3T922m2Vzbj4+M5SoqojezsgNGjpccDd+8CZ8/+\nGiI7d0pzZLm7q7dEAgIAe3vZSqc2eDBKSldsYRBRs2pqpGndG7ZEzp4FXFzUQyQwEHBwkLta0oSX\npIjIqO7flxaYahgip08DTk6Npz7p2VPuaqkhBgYRya6uTuoDaRgiubmAo2PjEHFykrta68XAICKT\nVF8P5Oc3noTRzq7pSRg59YnhMTCIyGwIARQUNJ76xNa28dQnAwcyRPSNgUFEZk0I4McfG4dIbW3j\nEPHwYIjogoFBRBapqUkYKyulEVkNQ8TTk5MwaouBQURW4+bNxvNn3brVeBLGIUM4CWNTGBhEZNVK\nSxuHSFERMGxY40kY27eXu1p5WW1gcD0MImpOWVnjSRgLC6VJFx+ehLFjR7mrNTyuh2F+ZRORjCor\nG0/CmJ8PeHk1noSxc2e5qzUMq21hmGHZRGRi7twBzpxRD5H//EfqSH94EsYuXeSuVncMDCIiPbp3\nT30SxpwcaT4td3f1EAkIALp1k7va1mFgEBEZWE2NNHNvwxA5cwbo37/xJIyOjnJX2zwGBhGRDO7f\nB/Ly1EPk1CmgT5/GU5/06iV3tRIGBhGRiairAy5dUg+RnBxp6veHQ6RvX+PXx8AgIjJh9fXSkrgP\n37XeuXPjEOnf37BTnzAwiIjMjBDAtWuNQ0ShaDx/lqur/kKEgUFEZAGEAH76qXGI3LvXOEQGDWpb\niDAwiIgsWFFR46lPysvVJ2EcMUK7SRgtMjCqqqqgVCoRHx+PiIiIRu8zMIjImpWUNA6RkpKmJ2Fs\n1+7X77PIwIiLi4O9vT28vb0ZGEaQmprKObn0iOdTv3g+tXP7duMQuX5dmurkQYjMn6/b306DzSIf\nHR0NJycn+Pn5qb2enp4Ob29veHp6IjExsdH3paSkwMfHB7179zZUafSQ1NRUuUuwKDyf+sXzqR1H\nR+Cxx4BVq4Ddu6U11q9fB956S+rz+OYb3ffRTvNH2iYqKgqxsbF44YUX1F5ftmwZtmzZAldXV0yc\nOBFz5sxBcnIycnJysGrVKqSlpaGqqgrnz59H586dMXnyZCi4xBYRUat16waMGyc9AOAf/9BtewYL\njNDQUBQUFKi9VlZWBgAICwsDAISHhyMrKwuRkZGIjIwEALz55psAgA8//BC9e/dmWBARmQphQFev\nXhW+vr6q5ykpKWL27Nmq55s2bRKvvfZaq7cLgA8++OCDjzY8dGGwFoYhCXZ4ExEZnVGXTg8KCkJe\nXp7q+blz5zBq1ChjlkBERG1k1MDo3r07AGmkVEFBAVJSUhAcHGzMEoiIqI0MFhhz5sxBSEgILl68\niAEDBuCDDz4AAGzYsAGLFy/G448/jpdeegm9WjHvr6YhuaSZm5sbhg0bhsDAQIwcORIAUFFRgenT\np2PgwIGYMWMGKisrZa7SdDU1XLyl8/fee+/B09MTPj4+yMjIkKNkk9XUuYyPj4eLiwsCAwMRGBiI\n5ORk1Xs8ly374YcfMH78eAwdOhRKpRKffPIJAD3/furUA2JkAQEBIi0tTRQUFIghQ4aIkpISuUsy\nO25ubuLWrVtqryUkJIiXX35Z3L17VyxZskS88847MlVn+tLT00VOTo7aYI7mzl9xcbEYMmSIuHbt\nmkhNTRWBgYFylW2SmjqX8fHxYt26dY0+y3OpWVFRkcjNzRVCCFFSUiLc3d1FeXm5Xn8/jXpJShcN\nh+S6urqqhuRS64mHBg0cP34c8+fPR8eOHREdHc3z2oLQ0FA4PrSkWnPnLysrC5MmTcLAgQMxbtw4\nCCFQUVEhR9kmqalzCTQ9qIXnUrO+ffsiICAAANCrVy8MHToU2dnZev39NJvAyM7OhpeXl+q5j48P\njh07JmNF5kmhUGDChAmYMWMGkpKSAKifWy8vLxw/flzOEs1Oc+cvKysL3t7eqs8NGTKE51YLiYmJ\nGDVqFBISElR/wI4fP85z2QqXL1/GuXPnMHLkSL3+fppNYJB+HD16FKdPn8batWuxcuVK3Lhxg8OU\nddSa88cbUVv24osv4urVq/jqq6+Qn5+PLVu2AGj6HPNcNq2iogLPPvss3n33XXTt2lWvv59mExgc\nkqsfzs7OAABvb29MmzYNBw4cQFBQEC5cuAAAuHDhAoKCguQs0ew0d/6Cg4Nx/vx51efy8vJ4bjXo\n06cPFAoFunfvjiVLlmDfvn0AeC61VVtbi5kzZyIyMhLTp08HoN/fT7MJDA7J1d2dO3dUTfySkhJ8\n9dVXmDRpEoKDg7Fjxw5UV1djx44dDOJWau78jRw5El999RUKCwuRmpoKGxsb2Nvby1ytaSsqKgIA\n3L9/H5988gkmT54MgOdSG0IIzJ8/H76+vli+fLnqdb3+fhqgs95gUlNThZeXlxg0aJDYuHGj3OWY\nnStXrgh/f3/h7+8vJkyYILZv3y6EEKK8vFxMmzZNDBgwQEyfPl1UVFTIXKnpmj17tnB2dhYdOnQQ\nLi4uYseOHS2evw0bNohBgwYJb29vkZ6eLmPlpufBuWzfvr1wcXER27dvF5GRkcLPz0+MGDFCrFix\nQm1EH89ly44cOSIUCoXw9/cXAQEBIiAgQCQnJ+v199Mk18MgIiLTYzaXpIiISF4MDCIi0goDg4iI\ntMLAICIirTAwiAwgLS0NmZmZcpdBpFcMDCID+Pbbb/Hdd9/JXQaRXnFYLZEWCgoKMHnyZIwfPx7f\nfvstwsLCkJiYiCeeeAIhISE4cOAAPDw8kJCQgE6dOmH06NGwtbVF7969kZiYiOrqavzxj39EWVkZ\nHBwckJaWJvchEbWaWS7RSiSHvLw8JCYmIjExEZMnT0ZmZiYUCgXy8/Nx8uRJ/Pvf/8aqVatw4MAB\nxMTEwN7eHitXrgQAKJVK7Ny5Ex4eHigvL5f5SIjahpekiLTUv39/PPbYY7CxscG4ceNUfRSzZ89G\nhw4dMHPmTOTk5KCmpgZCCLVJ38aOHYv58+dj586d6Nixo1yHQKQTBgaRlhwcHFRfd+jQAXfv3gWg\nPpPqg9k+H571880338TGjRtx/vx5+Pr6ora21ggVE+kXA4NIR5999hlqamqwb98+DB8+HB06dICr\nqytKSkpUn8nPz8ewYcOQkJCAjh07ori4WMaKidqGfRhEWmpurQAPDw+MGDECHh4eeOeddwAA4eHh\n+PjjjxEYGIjExES8++67uHTpEuzs7PDb3/4WLi4uxiydSC84SopIB+PHj8e6deswfPhwuUshMjhe\nkiIiIq2whUFERFphC4OIiLTCwCAiIq0wMIiISCsMDCIi0goDg4iItMLAICIirfw/2WuE7eTq6TEA\nAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## An illustration using matplotlib and scipy\n",
      "\n",
      "We define a function with a little more complex look"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def f(x):\n",
      "    return (x-3)*(x-5)*(x-7)+85\n",
      "\n",
      "x = linspace(0, 10, 200)\n",
      "y = f(x)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Choose a region to integrate over and take only a few points in that region"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "a, b = 1, 9\n",
      "xint = x[logical_and(x>=a, x<=b)][::30]\n",
      "yint = y[logical_and(x>=a, x<=b)][::30]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plot both the function and the area below it in the trapezoid approximation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(x, y, lw=2)\n",
      "axis([0, 10, 0, 140])\n",
      "fill_between(xint, 0, yint, facecolor='gray', alpha=0.4)\n",
      "text(0.5 * (a + b), 30,r\"$\\int_a^b f(x)dx$\", horizontalalignment='center', fontsize=20);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD3CAYAAADmBxSSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8jGf+//HX5BwycmgIRQgiCdJINIIV4tT6bg+01rYs\nug67lR4cqu23VRbdlm1V9bAr1G72u+qnVHXFWR0aWodMokVFCHKOEBIi50yS+f1xNyeCZExyT5LP\n8/GYx9y5Z+77/mQevOfOdV/3dWkMBoMBIYQQzYqF2gUIIYQwPQl3IYRohiTchRCiGZJwF0KIZkjC\nXQghmiEJdyGEaIbuGe7Tpk3Dzc0NX1/fO15bsWIFFhYWZGdnV6777LPP8PT0pFevXvz444+mr1YI\nIUSd3DPcp06dyp49e+5Yn5qayr59++jSpUvluszMTFatWsWBAwcICwtj1qxZpq9WCCFEndwz3IOD\ng3F2dr5j/WuvvcaHH35YY11UVBSjR4/G3d2doUOHYjAYyM3NNW21Qggh6sSqvhtERETQqVMnHnnk\nkRrrdTodPj4+lT97eXmh0+kYMWJEjfdpNBojSxVCiJatPgMK1OuCakFBAUuXLmXJkiV3HKy2g94t\nyA0GgzwMBhYtWqR6DebykM9CPoum8lkMHWoADHz4YeMet77qFe6XLl0iKSkJPz8/PDw8SEtLo1+/\nfly9epWgoCDOnj1b+d5z584RGBhY74KEEMJcRUXBoUPg6Agvvqh2NfdWr2YZX19frl69Wvmzh4cH\nJ06cwMXFhf79+/PGG2+QkpJCQkICFhYWaLVakxcshBBq+eAD5Tk0FNq0UbeW+7nnmfuECRMYNGgQ\n8fHxdO7cmX//+981Xq/e7OLm5kZoaCjDhw/npZde4tNPP22YipuRkJAQtUswG/JZVJHPooo5fRbn\nz8PWrWBrC7Nnq13N/WkMxjTmPMgBNRqj2o+EEEJNM2bAv/4Ff/4zrFnT+Mevb3ZKuAshxH1cvgwe\nHqDXK2fwnp6NX0N9s1OGHxBCiPtYuRJKSmDcOHWC3Rhy5i6EEPdw/Tp07Qr5+RAdDY8+qk4dcuYu\nhBAmtHKlEuz/8z/qBbsx5MxdCCHu4sYN6NIFcnPh6FEYOFC9WuTMXQghTOTTT5VgHzlS3WA3hpy5\nCyFELXJylLP2nBw4fBiCg9WtR87chRDCBD7/XAn2kBD1g90YcuYuhBC3uXVL6deenQ0HD8KwYWpX\nJGfuQgjxwFauVIJ98GDlzL0pkjN3IYSoJitLOWvPzVVGgBwyRO2KFHLmLoQQD+DDD5Vgf+wx8wl2\nY8iZuxBC/CojA7p3h8JC0OnAnKakkDN3IYQw0tKlSrCPHWtewW4MOXMXQgggOVkZFKy0FE6fhj59\n1K6oJjlzF0III/zlL8qQvs8/b37Bbgw5cxdCtHgnT0JAAFhZwblz0K2b2hXdSc7chRCint58EwwG\nePll8wx2Y8iZuxCiRdu7F0aPBkdHuHQJHnpI7YpqJ2fuQghRR2Vl8MYbyvI775hvsBtDwl0I0WKt\nWwe//KKM/vjqq2pXY1oS7kKIFik3VzlbB3j/fbCzU7ceU7tnuE+bNg03Nzd8fX0r173xxhv4+PgQ\nEBDAnDlzKCwsrHzts88+w9PTk169evHjjz82XNVCCPGAli5V7kjt3x8mTFC7GtO7Z7hPnTqVPXv2\n1Fj32GOPERsbS0xMDPn5+WzYsAGAzMxMVq1axYEDBwgLC2PWrFkNV7UQQjyAixfh44+V5c8+A4tm\n2IZxz18pODgYZ2fnGutGjRqFhYUFFhYWPP744xw6dAiAqKgoRo8ejbu7O0OHDsVgMJCbm9twlQsh\nhJFeew1KSuCFFyAoSO1qGobVg2y8du1aZsyYAYBOp8PHx6fyNS8vL3Q6HSNGjLhju8WLF1cuh4SE\nENJUB0wWQjQ5e/fC9u2g1cKyZWpXc3eRkZFERkYavb3R4f7uu++i1WoZP348QK39LzUaTa3bVg93\nIYRoLCUlMGeOsrxwIXTooG4993L7ie+SJUvqtb1RLU3/93//x969e1m/fn3luqCgIM6ePVv587lz\n5whs6sOqCSGalRUrlOEFPD1h9my1q2lY9Q73PXv2sHz5crZt24Zdtb5D/fv3Z+/evaSkpBAZGYmF\nhQVardakxQohhLESEuDdd5XlVavAxkbdehraPZtlJkyYwKFDh7h+/TqdO3dmyZIlLFu2jJKSEkaO\nHAnAwIEDWbVqFW5uboSGhjJ8+HBsbGxYs2ZNo/wCQghxPxXjxhQVwcSJ8Gt8NWsytowQotn7+mt4\n7jlwclKaZdzc1K6o/mRsGSGEqCYnp+oi6t/+1jSD3RgS7kKIZu2114rJyICBA+FPf1K7msYjzTJC\niGbrv//N4dlnHbG0LOPnny2pNpJKkyPNMkIIAVy8eKXyTH3UqKgmHezGkHAXQjQ7SUlJTJ9+haws\nR3r2zGfEiGi1S2p0Eu5CiGbl/PnzrFgRzeHDfbG0NLBgQQKWluVql9XoJNyFEM3GqVOn2Lp1P//9\n75MATJ2aQc+eBSpXpY4HGjhMCCHMxfHjx4mKiuKHH54jPd0eT88Cpk+/onZZqpFwF0I0aQaDgUOH\nDnHmzBny8x9n58722NiU8957iVhbGygtVbtCdUi4CyGarLKyMr777jsSExPp2nUQkyb1AGDWrDS6\ndy9SuTp1SbgLIZokvV7Pzp07yczMJCDgUebM6UFOjhUDB+bw3HPX1C5PdRLuQogmp7i4mK1bt5Kf\nn09AQADr1rUnKqoNTk56Fi1K4i5TSbQoEu5CiCYlPz+fLVu2YGFhgZ+fHydPtmbVqo4ALFqUjKtr\nC21kv42EuxCiycjJyWHz5s04ODjg6enJjRtWvP12N8rKNEyZcoXg4By1SzQbEu5CiCbh+vXrbN68\nGTc3N7p06UJZGSxc2JVr12zw88vjpZfS1S7RrEi4CyHMXkZGBlu2bKFLly48/PDDAPzrXx04ftwR\nJyc9S5cmYCVpVoN8HEIIs5acnExERAQ9e/akbdu2AERGOvLFFw+j0Rj461+TcHPTq1yl+ZFwF0KY\nrfPnz7N792769OmDs7MzAJcu2fGXv3gA8Mor6QwceEvNEs2WhLsQwiydOnWK77//nr59+6LVagHI\nybFk3rzuFBRY8vjj2UyZclXlKs2XhLsQwuxUjBPTr18/WrVqBUBpKbz9djfS0uzw8ipg4ULpz34v\nEu5CCLNRMU7ML7/8QmBgILa2tr+uh6VLu6DTtcHZWc+KFRexs5MZ3e5Fwl0IYRaqjxPTv39/rKp1\nfwkPb8+2ba7Y2pazcuVF2reXC6j3c8/x3KdNm4abmxu+1eanys3NZcyYMbi7uzN27Fjy8vIqX/vs\ns8/w9PSkV69e/Pjjjw1XtRCiWdHr9Wzfvp3U1FQeffTRGsG+c6cLYWEd0WgMLF2aQJ8+LXN89vq6\nZ7hPnTqVPXv21FgXFhaGu7s7Fy5coFOnTqxevRqAzMxMVq1axYEDBwgLC2PWrFkNV7UQotkoLi7m\n22+/JTs7m4CAACwtLStfO35cy7vvdgXg9ddTGTpU7kCtq3uGe3BwcGX3owo6nY7p06dja2vLtGnT\niIqKAiAqKorRo0fj7u7O0KFDMRgM5ObmNlzlQogmLz8/n02bNqHX6/Hz80NT7QrpyZOtmTevB2Vl\nGiZNuiIjPdZTvdvco6Oj8fb2BsDb2xudTgco4e7j41P5Pi8vL3Q6HSNGjLhjH4sXL65cDgkJISQk\npL5lCCGauNvHianu3Dl7Zs/2pLjYgqefvs6sWS1vaIHIyEgiIyON3r7e4W4w1P0KteYu/ZSqh7sQ\nouW5fZyY6hIT7XjlFU/y8y0ZOTKbd95JxqIFzvZ8+4nvkiVL6rV9vT+ywMBA4uLiAIiLiyMwMBCA\noKAgzp49W/m+c+fOVb4mhBAVMjIy2LhxIx07drwj2C9dsmPmzJ7cvGnNb36Tw1//mkS1JnhRD/UO\n96CgIMLDwyksLCQ8PJwBAwYA0L9/f/bu3UtKSgqRkZFYWFhU3lUmhBCgjBOzefNmevToQceOHWu8\nFh9vz4sv9iQry5r+/W/xwQeXsLaWvuzGumezzIQJEzh06BBZWVl07tyZd999l9DQUCZNmoSXlxcB\nAQF88MEHALi5uREaGsrw4cOxsbFhzZo1jfILCCGahtrGialw7pw9L7/ck5wcKwYNymH58kvY2kqw\nPwiNoT6N6KY4oEZTr3Z7IUTTd/r0aQ4ePFhjnJgKJ044MG9ed/LyrBgy5CZ/+1sCNjamy4jS0lKO\nHj3K3LlzTbZPNdQ3O+UOVSFEg6ptnJgKBw44sXChByUlFowYcYP33kuUphgTkXAXQjQIg8HA4cOH\nOX36dI1xYips2tSWjz7qjMGgYfz4TF5/PVUunpqQhLsQwuTuNU5MaSl8/nkn/t//cwPgpZfSmTr1\niozwaGIS7kIIk9Lr9ezcuZOrV6/y6KOP1hhOIDfXkvnzPTh2zBFLSwPvvJPM009nqVht8yXhLoQw\nmeLiYrZu3Up+fj79+vWrcSNjUpIt8+b1IDnZDicnPR9+mEBAQN499iYehIS7EMIk8vPz2bJlCxqN\nBj8/vxqv7d3rzPvvd6GgwBJPzwJWrLjEww+XqFRpyyDhLoR4YBXjxGi1Wnr06FG5vqhIw8cfd+bb\nb5WJrUeNymbhwmRatSpXq9QWQ8JdCPFAsrKy2Lx5M23btqVr166V6y9csGfhwq5cvNgKG5ty5s1L\n5dlnr8uF00Yi4S6EMFpGRgZbtmzB3d29cjiB0lL4z3/as3ZtB0pLLXB3L2LZsgS8vApVrrZlkXAX\nQhglOTmZrVu30rNnT9q1awfAxYt2vPtuV86ebQ3A+PGZvPpqujTDqEDCXQhRb/Hx8ezatatynJj8\nfAu++OJhNm5sR1mZhvbti1m4MJmgIJmwRy0S7kKIeqk+Tkzr1lr27nXmk086ce2aDRYWBn7/+0xe\neikdBwc5W1eThLswWnk5pKbC+fOQlARXryqPK1eU5+xsKCys+dDrwdISrKyqHjY2Blq1KkOrLcfZ\n2QIXF0ucnTW4uUH79tChg/L88MPQuTPcdhe7aEQ6nY5jx47Rr18/YmPd+OyzjsTFKU0wvXvn89Zb\nKfj4yATW5kDCXdTJjRsQHa08Tp1SAv3CBSWw66u0VHlU0VDXf4oaDXTsCB4eVY9u3aB7d/DxAReX\n+tcj7q/6ODH29sN4660uHD3qCICrawkzZ17m6aezWuSMSeZKwl3UKjMTIiPh4EH4/nuIj6/9fW5u\n4OUFPXooZ9dublWPhx6CVq3A3h70+hyysy+Tnp5IQkIyJSXltG7tSJs2LrRp0xa93p68PEvy8iy5\ndQuuXSvj6lUDmZmWZGXZcPOmHTdutCInR0tamgVpafDDD7XX4+MDvXopj4plNzekC56RysvL2bv3\nO/bsKebIkSmcOKGEeuvWZUyZcoWJEzOxt5cmGHMj4S4AMBggNhYiIpRHdHTN1+3sICAAAgOVZ29v\nJdQdHWvfX35+Punp6SQkJJOQkEBhYSFt2rTB2dmZgABf7O3tb9uiuE515uUVk5CgJykJUlOtyciw\nJzOzFZmZTmRmunD1qjVXrypfTNU5O0Pv3uDrC336VD3fNmeEuE1Ojp7588+wffujpKa6AkqoP/dc\nJhMnXsXJqUzlCsXdyGQdLVxSEqxfrzzOn69ab2cHgwfD8OHKIyAArK3vvp+ioiIuX75MSkoKCQkJ\n5OTk4ODggLOzM23btsXBwaFBf4+ysjJyc/NJTCwlPt6KlJTWpKY6cPmyE1evulBYaFfrdh071gx7\nX1/lbP+O754WpLzcgE6nZ926ctatsyA/3wYAJyc9EyZk8vvfX0OrbTqh3lIn65Bwb4FKSmDLFli9\nGg4frlrv6gpPPQVjxsCoUUqTyt33UUJGRgapqalcunSJrKwsHBwccHR0xNXVFUdHxxqDRqmpqKiY\nlBQ9cXEWXLxoT3KylrQ0JzIyXNDr7/zGsrAw0KOH5o7Q795duQDcVJSVlVFcXFzro6ioiMLCQgoK\nCigsLKSwsJCEBCt0uk7odD25etW1cj99+uQxfvw1Ro680SSnvpNwbyQS7uq5cgVWrYIvvlB6s4By\nhvrMMzB5MowceffwKisr48qVK6SlpXHp0iWuXLlC69at0Wq1tG3bFmdnZ7MJ87oqLTVw6VIZZ89a\nEh9vQ1KSltRURzIznSgvv/PKoK2tAR8f8PXV1Aj9jh0btj2/pKTkriF9e0BXLBcVFaHX67GxscHS\n0hIrKyssLS0rH8r46rYkJrpx6lQHoqLak5pa9deVk5Oexx67wZNPZtGrV9Pu/dJSw70JnYcIYyUl\nwfLl8K9/QfGvTdt9+sDLL8PEidCmzZ3blJeXc+3aNdLS0khISCAtLQ07Ozu0Wi2urq706NGjxjjd\nTZGVlQYvLyu8vAD0QDaQTXGxhkuXrIiNtSA+3oaEhNakprYhO1vLyZNw8mTN/Tg6ltOrVzl+fpb4\n+lad8Vdvzy8vL7/vWXRFOBcUFFSuKyoqQqPRYGlpibW1dY1wrni2sbHB2toarVaLi4sLtra2WFtb\n15ggAyAnx5K4uFbExrbmxAktp045UFxc9SWm1ZYyeHAOo0bdYNCgnCb1V4q4k5y5N2OpqbB4Maxb\nV9X18JlnYM4cCA6uebZpMBjIzs4mPT2dxMREkpOTsbS0xMHBAVdXV1xdXe8Ii5YmL8+CS5fsOHvW\nivh4ay5dakVychvy82tvz3d0LMTF5RaOjjk4Od3ExSWPtm3zadu2ABeXIrTaMqytLSuD28bGpjKo\nqz9b1KN/YWkpXL9uTWamDenptiQm2pGQYMfFi/akpd1Zp6dnAYGBuQwZcpO+ffOaZaDLmbtoNm7c\ngGXL4LPPlDN1S0uYNAneflvpFljh5s2bXL58maSkJBITEykvL688+3v00UfvmPOypXNwKMfPr4Cq\nocqzMBggK8uKS5fsOX/elvPnrUlIaEVKipacHHtycuwBt1r3Z2lpwMmpFGdnPS4upTg4lNGqVRn2\n9uWVD2vrO7sYlpRYUFBgQUGB0nW0oMCCmzetyMy04fp1a8rLa28jsrUtx8urgN6983nkkXwefTQX\nZ+fSWt8rmj6jw33t2rX8+9//pri4mODgYD755BNyc3OZNGkSP//8MwEBAaxfv77Be0mIKmVlSnv6\nO+8oAQ/w3HPw3ntKP/S8vDzOn08nOTmZxMRECgsLK8Pcz8+vlu6J4n40GnB1LcXVNbfGOCplZXDt\nmjUZGbZkZNiQkWHDlSs2vy7bkpVlRV6eFVlZ1mRl3aMbUr3rMeDqWoKbm5727Uvw8CjEw6OIbt2K\n8PAobJZn5qJ2RjXLZGdn069fP86cOYO9vT1PPvkks2fP5tSpU6SmpvLRRx8xb948unbtyuuvv17z\ngNIs0yCio+GllyAmRvl52DB4990iOnRIJzU1lYsXL5Kbm4tWq8XJyalRuieKeysp0XDzphXZ2Vbc\nvGlFbq4VhYUW1R6W6PWaOy7WWlkZaN26rPLRqlU5bdqU4uamx9VVj7W1/P+qTppl6sHe3h6DwUBO\nTg4ABQUFODk5odPpWLBgAba2tkybNo1ly5YZs3tRD/n5SnPL3/+u3IjUoUMpL798gY4ddeh0Vd0T\nu3XrRps2bZpcj5bmzMbGQLt2etq106tdimiGjA73sLAwunbtiq2tLbNmzSIoKIjo6Gi8vb0B8Pb2\nRqfTmbRYUdOPP8LkyaUkJVlhaVlOcPAJxo49Tbt2rXB1fZg+fXpLmAvRQhkV7teuXSM0NJSzZ8/i\n7OzM+PHj2bFjR53/ZFi8eHHlckhICCEhIcaU0WIVFyvt6h9/bMBgsKJz5xv87//GERhohaVlX7XL\nE0KYQGRkJJG3j6NRD0aFu06nY8CAAZUT4Y4fP54ffviBwMBA4uLi8Pf3Jy4ujsDAwFq3rx7uon4S\nEuD3v4cTJ5Q7KZ9/PoFZs3KwtpaeLUI0J7ef+C5ZsqRe2xs1QGdwcDAxMTFkZ2dTXFzM7t27eeyx\nxwgKCiI8PJzCwkLCw8MZMGCAMbsXd7FlC/j7K8Hetm0+ixcfYN68m3IBTQhxB6PCvU2bNixYsIBn\nnnmGwYMH4+fnx7BhwwgNDSUlJQUvLy/S09OZOXOmqettkcrK4I034He/g1u3IDj4Om+99TW//a0M\nXi6EqJ3coWrmbt6ECRNgzx5l3Je33rqOi8t6Bg4cgPW9hmkUQgAttyukzJtixuLjYcAAJdhdXWHr\n1lzatfuKRx7xlWAXQtyThLuZOnoUBg5Uxlh/5BE4dqyU7OytPPzwwzg5OaldnhDCzEm4m6Ft22DE\nCGWC6aeegiNHIDHxe0pKSujatava5QkhmgAJdzOzdq0ycmNREfzpT/Dtt5CScpa4uDj69OmjdnlC\niCZCwt2MrFgBf/4zlJfDokWwZg3cuHGNffv24efn1+KH3BVC1J2khZlYtgzmz1eWV62C0FAoLi4m\nIiICDw8PGeRLCFEvEu4qMxjg3XeVSTU0GmW2pKlTlckzvvvuO2xsbOjYsaPaZQohmhhpllHZkiVK\nsFtYKDMmTZ2qrP/pp59ISUnBx8dH1fqEEE2ThLuKVqxQwt3CAjZsUGZLArh8+TI//PADfn5+9Zpi\nTQghKkhyqGTtWqiYxyQ8XJkxCSA/P5+IiAi8vLxkZiQhhNEk3FWwcSO8+KKy/Pnn8MILynJ5eTm7\ndu3CycmJdu3aqVegEKLJk3BvZAcPwpQpyoXUpUvhlVeqXjt27BhZWVl4enqqV6AQolmQcG9EZ87A\ns8+CXg9z5yrT41VITEwkOjoaPz8/mT1JCPHAJNwbyeXL8NvfQk4OjBsHH31U9VpOTg47d+6kT58+\n2NjYqFekEKLZkHBvBLm58MQTkJoKgwbBl18qPWQAysrK2LFjB+3bt8fZ2VndQoUQzYaEewMrL1cu\nmJ48CZ6eEBEB1TvBREZGUlhYiIeHh3pFCiGaHQn3BvbXv8J//wuOjrBjhzIue4Vz585x5swZfH19\n1StQCNEsSbg3oP/+t+ru040boWfPqteysrLYu3evDAgmhGgQEu4N5MwZmDxZWf7b32D06KrXSkpK\niIiIoGvXrmi1WnUKFEI0axLuDeDWLaXLY34+TJxYdSdqhX379mFpaUmnTp3UKVAI0exJuJuYwaCM\nyX7hAvj6KsMMVO+2fvLkSRITE2VAMCFEg5JwN7HVq2HTJnBwgM2boVWrqtcyMjI4dOgQfn5+WFpa\nqlekEKLZMzrc8/PzeeGFF+jZsye9evUiKiqK3NxcxowZg7u7O2PHjiUvL8+UtZq9n36COXOU5X/+\nE7y8ql4rLCxk27Zt9OjRg1bVE18IIRqA0eG+aNEi3N3dOX36NKdPn8bb25uwsDDc3d25cOECnTp1\nYvXq1aas1azdugXjx0NJiTKLUsUoj6BMvLFr1y4cHBxo3769ekUKIVoMo8N9//79zJ8/Hzs7O6ys\nrHB0dESn0zF9+nRsbW2ZNm0aUVFRpqzVrM2aBQkJ0LcvfPxxzdeioqLIzMzEq/qpvBBCNCCjwj0t\nLY2ioiJCQ0MJCgrigw8+oLCwkOjoaLy9vQHw9vZGp9OZtFhztXkz/Oc/YGenTLphZ1f1WnJyMseP\nH5cBwYQQjcqou2eKioqIj49n+fLljBw5khdffJGvv/4ag8FQp+0XL15cuRwSEkJISIgxZZiFtLSq\nsdlXrIDqnWByc3PZsWMHvXr1wtbWVp0ChRBNUmRkJJGRkUZvrzHUNZFv4+PjQ1xcHAC7d+9m3bp1\nlJSUsGDBAvz9/Tlx4gTLli3jm2++qXlAjabOXwLmrrwcRo1Sxmh/4gnYvr2q22NZWRmbNm3CysqK\n7t27q1uoEC1YaWkpR48eZe7cuWqX8kDqm51Gt7l7enoSFRVFeXk5O3fuZOTIkQQFBREeHk5hYSHh\n4eEMGDDA2N03CZ99pgR727bwr3/V7M/+448/kp+fL8EuhFCF0eH+0UcfMXv2bAICArCzs+P5558n\nNDSUlJQUvLy8SE9PZ+bMmaas1axcvAjz5yvL//wnuLlVvRYfH8/Jkyd55JFH1ClOCNHiGT1iVc+e\nPTl+/Pgd6yMiIh6ooKagvBymT4fCQvjDH+Dpp6tey87OZvfu3fj6+sqAYEII1cgdqkZYtQoOH1bO\n1j/9tGq9Xq9n+/btuLu74+joqF6BQogWT8K9nhIT4a23lOVVq+Chh6pe279/P+Xl5bi7u6tTnBBC\n/ErCvR4qBgXLz1fuQH322arXfvnlFy5evEjv3r3VK1AIIX4l4V4PGzbA/v3g4gKff161/urVqxw4\ncEAGBBNCmA0J9zq6cQNee01ZXr5c6f4Iyg1dERER9OjRg9atW6tXoBBCVCPhXkdvvw2ZmRAcDFOn\nKusMBgN79uzB3t6eDh06qFugEEJUI+FeB8eOwZo1YG2tjNdecbNSTEwMly9flok3hBBmR8L9PvT6\nqrFj3ngDevVSllNTUzly5IgMCNaI1q1bR3BwMGfOnFG7FCHMnoT7fYSFwS+/QLdusGCBsi4vL4/t\n27fj4+ODXfUhIEWDGjduHPb29tIjSYg6kHC/h2vXYNEiZXnlSrC3h/Lycnbs2IGLiwuurq7qFtjC\nxMTE4O/vL38pCVEHEu73sGAB3LwJjz8OTz2lrDty5Ag5OTl4enqqW1wLFBUVhVar5fDhw/ztb3/j\n4sWLapckhNmScL+Ln3+GtWvBygo++US5iHrp0iVOnDiBn5+f2uU1e4cOHeKZZ55h+vTpJCcnA0q4\njxkzhiFDhjBo0CBWrVqlcpVCmC8J91oYDPDqq8rzrFng7Q03b95k586d+Pr6Ym1trXaJzdrZs2d5\n8803WbJkCYWFhaxYsYIrV65gMBjw9fUFlBvHCgoKVK5UCPMl4V6LTZvgyBFo1w7+8hdlQLBt27bR\nqVMnnJyc1C6v2fv888/p378/vX7tmtShQwfOnTtHnz59Kt9z/PhxAgMD1SpRCLMnY9Lepri4amCw\n998HR0fYty8SvV5Ply5d1C2uBYiNjSUmJoa3334bKysrNmzYAMCFCxcqv1hTUlJISkri/fffV7NU\nIcyahPs1Tm1gAAAUnUlEQVRt/v53SE6GPn2UO1FjY2OJi4sjKChI7dJahL179wIwdOjQGus9PT1p\n164dERERJCQksGbNGumGKsQ9SLhXk50N772nLH/4IWRnX2P//v307dtXJt5oJAcOHMDDw4OHqo+l\n/KtJkyapUJEQTZO0uVfz/vtK18cRI2DYsGIiIiLw8PDAwcFB7dJahOTkZDIzM+nbt6/apQjR5Em4\n/yoxUWmSAfjwQwP79n2HjY0NHTt2VLewFiQmJgagxoVTIYRxJNx/9c47UFICkyaBwfATKSkpMiBY\nIztx4gSAfO5CmICEO3DqFHz1FdjYwCuvXOHw4cP07dsXCwv5eBrTiRMnsLGxoVu3bmqXIkSTJ+kF\nLFyoPM+Yoeenn/4rA4KpICkpiezsbLp16yazWQlhAkaHe1lZGf7+/jz166Arubm5jBkzBnd3d8aO\nHUteXp7JimxIx4/D9u3QqpWBvn134+TkRNuKaZZEozl58iQAPXv2VLkSIZoHo8P9008/pVevXpUj\n9IWFheHu7s6FCxfo1KkTq1evNlmRDaliGN9x49IoK7ssA4Kp5KeffgIk3IUwFaPCPS0tjV27djFj\nxgwMBgMAOp2O6dOnY2try7Rp04iKijJpoQ3h++/hwAHQasvw9t4hE2+o6JdffgGgR48eKlei/FVq\nrNLSUhNWIoTxjAr3uXPnsnz58hoXHKOjo/H29gbA29sbnU5nmgobiMGg9JABGDIkmv79PbGxsVG3\nqBbqxo0bpKWlodFo6N69u6q1xMTEsHXrVqO3X716deUolkKoqd63Xe7YsYN27drh7+9PZGRk5fqK\nM/i6WLx4ceVySEgIISEh9S3jge3dq8yNqtUW8fzzV3F27tToNQjF6dOnAXB2dm6UgdlSU1MJCwuj\nbdu26PV63nzzTQDOnDnD7t27WVhxhd0IkydPZs6cOaxcubLOv8vKlSvZu3cvWVlZrF69mn79+hl9\nfNF8REZG1sjY+qp3uB89epRt27axa9cuioqKuHXrFpMnTyYwMJC4uDj8/f2Ji4u754h91cNdDQYD\nLFmiLD/++Cl8fCTY1dSYTTJ6vZ5XXnmFGTNm8Msvv7Br1y5mz54NwPLly1mzZs0D7d/R0ZHf/e53\nzJs3jy+++KJOPX/mzp1Lx44d+fTTTyuHNBbi9hPfJRWhVUf1bpZZunQpqampJCYmsnHjRoYPH86X\nX35JUFAQ4eHhFBYWEh4ezoABA+q760Zz4IDSS6Z160Jeekl6g6qtItwb42L2sWPHuHz5MgEBAYwZ\nM4awsDBsbW356quvGDx4sEm6wD7xxBNYWVlx6NChOm9z8uRJevXqJU2DwmQeONkqLkCGhoaSkpKC\nl5cX6enpzJw584GLayh/+YsegIkTM2jTRsJdTWVlZZw9exZonHA/ceIETk5OdOzYkd69e+Pr60tx\ncTHr16/nd7/7ncmO8/LLL7Nly5Y6v//nn38mICDAZMcX4oGGOhw6dGjl0KxarZaIiAiTFNWQDh2C\nY8esad26mD/8IUftclq8xMREioqK0Gg0jRLusbGx9O7du8a6mJgY2rdvj7Ozs8mO0717d2JiYkhL\nS6NTp3s3+6WlpXH9+nUJd2FSLW4c23ffVZ5HjTqLg0O5usUI4uLiALCysmrQYQeWLl3KlStXOHXq\nFF27dmXWrFm4u7vz+uuvc/To0XvOi5uQkMCOHTsoKSkhLy+P+fPn8+WXX5KTk0NWVhavvvoq7du3\nr7FN69atcXFx4dChQ/zhD3+o8dq5c+c4ePAger2enJwcvLy8sLS0vKMGY44rRIUWFe4//ggHD4KD\nQxmjRsUBXmqX1OJVNMl4eHg06Jj58+fPJz09nbFjx/Lyyy/XuFB19uxZnn766Vq3y8jIICIigrlz\n5wLw1ltvMXnyZObNm4dWq2Xq1KkEBgYyduzYO7bt0qULly9frrHu+PHjLFq0iPXr19O2bVuSkpKY\nMGECvXv3rtHe/yDHFQJa2NgyS5cqz5MmZdOqVYm6xQigKty9vBr+i/b8+fPAnXfBZmdno9Vqa93m\n66+/rnH9SK/XY2dnR//+/XFxcWHatGmMHDmy1m3d3d3JyMio/Pny5cu88847zJ07t3KIi65du9Kq\nVas7mmQe5LhCQAsK99OnYfdusLeHyZNvqF2OQLmYevHiRaBxhvmNj4/HwcGBhx9+uMb6e4X7+PHj\nsbe3r/w5Li6usieYm5sbf/7zn+86mUuXLl24cuVK5c+ff/45paWlDB8+vHJdQkICt27duiPcH+S4\nQkALCvcPP1SeZ8wAZ2fjby8XppOUlERJSQkajabRwr22sWs0Gg35+fm1blP9iyApKYlr167x6KOP\n1ul4ZWVllJcr13UMBgMxMTEMHDiwRnfHEydOYGFhccfsUw9yXCGghYR7UhJs3AiWlvDaa2pXIyrE\nx8cDysXUiqErGvp4tTX/ODs7k5SUdN/tY2JisLa25pFHHqlcl5aWdtf3JycnV84Fm5SUxM2bN+/4\ncomJicHHxwd7e3vS09NNclwhoIWE+8cfQ1kZTJgAXbuqXY2ocOHCBUC5M7WhJyC/efMmV69erbW7\npaurKykpKXesLykpYe3atZVNR0ePHsXDwwNbW1sACgoK+Prrr+96zOrh3rZtW6ytrencuXPl60VF\nRfz000/4+/sD8NVXX5nkuEJAC+gtc/06/POfyvKvQ4gIM1ERXo0xZ2rFxdTawt3X15dTp07dsf7E\niRN88cUX9OzZk9LSUq5cuVL5JaTX6/nnP//JxIkT73rMlJQURo8eDYCDgwOBgYGVXyKlpaWsWLEC\ngIceeogrV67QoUMHkxxXCGgB4f73v0NhIfz2tyDDdpiXinC//aaihnD+/Hm0Wm2tbe4DBw5k27Zt\nd6z39fVl9OjR6HQ6bGxsWLduHZ988glLly5Fq9UyevTou/Yzv3XrFjdu3GDQoEGV6xYuXMjGjRv5\n8MMPKSsr48UXX+Q3v/kN69atIy8vjylTpjzwcYWo0KzDvaBACXeA//1fdWsRNeXm5nLt2jU0Gk2j\nhPu5c+cIDAysdV5cf39/LCwsuHz5co0LmQ4ODvz1r3+t8d7XX3+9Tsc7f/48PXv2rLE/V1dXXnnl\nlRrvq21U1Ac5rhAVmnWb+5dfQlYW9O8PwcFqVyOqu3TpEgBt2rShawNdCPn222+ZNWsWoPSn/+1v\nf1vr+2xsbJg+fTqffPKJSY5bXl7O3//+d1588UWT7E8IYzTbcC8vh4r/q3PngkywZF4SEhIA7ugC\naEo7d+7EycmJM2fO4OLiUjkOUm3Gjx9PfHw8P/zwwwMfd/PmzVhbWzNkyJAH3pcQxmq24b53L5w7\nB506wbhxalcjblcR7hU9RRrClClTsLW15cCBA3c0c9zOysqKjz76iNWrV1NUVGT0Ma9du8a3337L\ne++9Z/Q+hDCFZtvmvnKl8vzKK2BtrW4t4k4V3SAb8sy9+qilddGjRw/efvttNm3axAsvvGDUMTds\n2MDy5cvlgqdQXbMM9zNnYN8+aNUK/vxntasRtblw4QL29vaNcvNSffTp0+eBumZWzOokhNqaZbNM\nRVv7H/8IJhyiW5hIRkYGubm59OnTp07T0Akh6q/Zhfu1a7B+vbIsJ1HmqWIkSJkIWoiG0+zCfe1a\nKC6GJ5+EWu5XEWagItzNeZ5dIZq6ZhXupaUQFqYsv/qqurWIuzt16hROTk6NcvOSEC1Vswr37dsh\nLQ08PUHmMTBPBQUFnDlzhqCgILVLEaJZa1bhXjHUwMsvQy13mQszEB0dTVlZGcFyy7AQDarZRGBc\nnDI/aqtWYGQXZdEAvvjiCyZMmEBpaSmgDAnQsWNHRo0apXJlQjRvRoV7amoqw4YNo3fv3oSEhLBh\nwwZAGQxqzJgxuLu7M3bsWPLy8kxa7L384x/K8+TJ4OTUaIcV93H06FE0Gg0ajYa0tDSOHz/On/70\np1oH8BJCmI5R/8Osra1ZuXIlsbGxfPPNNyxYsIDc3FzCwsJwd3fnwoULdOrUidWrV5u63lrdugX/\n+Y+y/PLLjXJIUUfDhg2jc+fOnDt3jnnz5uHp6XnXAbyEEKZj1B2q7du3r7y92tXVld69exMdHY1O\np2PBggXY2toybdo0li1bZtJi7+bLLyEvD4YMkTHbzc2zzz7L1atXmTt3Lv369WP+/Plo7jKKm8Fg\nYMOGDTg6OpKVlUVqaip//OMf6dSpUyNXLUTT98B/G1+8eJHY2Fj69+9PdHR05e3k3t7e6HS6By7w\nfgwGqPgDQc7azY9Wq+XNN9/ku+++Y9myZWi12ru+NywsDAsLC5588knGjBnDwYMHJdiFMNIDjS2T\nm5vLc889x8qVK3FwcMBgMNRpu8WLF1cu1zZZQX0cP66MJdOuHYwda/RuhMrS09P56quv+O677wDl\npCEgIEDlqoRQT2RkJJGRkUZvb3S46/V6xo0bx+TJkxkzZgwAgYGBxMXF4e/vT1xcHIGBgbVuWz3c\nH9QXXyjPU6eCjY3JdisaWXR0NL6+vtjb2wOg0+kIDAwkNzf3nmf7QjRXt5/4LlmypF7bG9UsYzAY\nmD59On369GHOnDmV64OCgggPD6ewsJDw8PAGv7385k3YtElZnjGjQQ8lGljbtm1p164doNzo9P33\n3/PII4+wf/9+lSsTomkyKtyPHDnC+vXrOXjwIP7+/vj7+7Nnzx5CQ0NJSUnBy8uL9PR0Zs6caep6\na1i/Xpn8esQI6NGjQQ8lGtiAAQPo0KED+/fv58KFCzz77LMcOHCALl26qF2aEE2SUc0ygwcPpry8\nvNbXIiIiHqigujIYqppkZMz2ps/S0rLGnKN+fn4qViNE09dk7ySJioJffoG2beVCqhBC3K7JhnvF\nWfsf/ygXUoUQ4nZNMtxzcmDjRmX5T39StxYhhDBHTTLcN29WLqQOHaoM7yuEEKKmJhnu//d/yvPU\nqaqWIYQQZqvJhXt8PBw5Aq1bw7hxalcjhBDmqcmFe8Xoj+PHg4ODurUIIYS5alLhXlYG69Ypy3/8\no6qlCCGEWWtS4X7woDJHarduILO0CSHE3TWpcK+4kPrCCzJHqhBC3EuTicicHPj2W2V5yhR1axFC\nCHPXZMJ90yYoKoJhw6BrV7WrEUII89Zkwl36tgshRN01iXA/fx6OHVO6Pj77rNrVCCGE+WsS4V7R\nt/33v1duXhJCCHFvZh/u0rddCCHqz+zD/dAhSE9X+rYPHqx2NUII0TSYfbhv2KA8T5wIGo26tQgh\nRFNh1uFeXAzffKMsT5yobi1CCNGUmHW4796t3LzUty/4+KhdjRBCNB1mHe5ffaU8T5igbh1CCNHU\nmG245+bCtm3K8vPPq1uLEEI0NSYP98OHD+Pj44Onpyeff/650fvZulUZbiA4GNzdTVigGTlx4oTa\nJZgN+SyqyGdRRT4L45k83GfPns2aNWvYv38///jHP7h+/bpR+2kJTTLyD7eKfBZV5LOoIp+F8Uwa\n7jk5OQAMGTKELl268NhjjxEVFVXv/Vy7Bt99B1ZWyoxLQggh6sfKlDuLjo7G29u78udevXpx/Phx\nnnjiiXrtZ/Nm5c7U//kfcHU1ZYUKjUZDdnY2P/30k+l3Xg8ZGRmq12Au5LOoIp9FFVN8FuXl5VhZ\nmTTqmgRVfmNNHe9G2r27+d+4tH37drVLMBvyWVSRz6KKqT6LWbNmmWQ/TYVJwz0wMJA33nij8ufY\n2FhGjx5d4z0Gg8GUhxRCCFELk7a5Ozo6AkqPmaSkJPbt20dQUJApDyGEEKIOTN4s88knn/Diiy+i\n1+uZNWsWrg3RaC6EEOKeTN4VcujQocTFxXHx4sUabVym6v/eHKSmpjJs2DB69+5NSEgIGypGR2uh\nysrK8Pf356mnnlK7FFXl5+fzwgsv0LNnz8rOCC3V2rVrGTRoEP369WPOnDlql9Oopk2bhpubG76+\nvpXrcnNzGTNmDO7u7owdO5a8vLz77qfR7lA1Vf/35sDa2pqVK1cSGxvLN998w4IFC8jNzVW7LNV8\n+umn9OrVq84X2purRYsW4e7uzunTpzl9+jQ+LXRApezsbJYuXcq+ffuIjo4mPj6evXv3ql1Wo5k6\ndSp79uypsS4sLAx3d3cuXLhAp06dWL169X330yjhbqr+781F+/bt6du3LwCurq707t2bmJgYlatS\nR1paGrt27WLGjBkt/mL7/v37mT9/PnZ2dlhZWVVew2pp7O3tMRgM5OTkUFhYSEFBAc7OzmqX1WiC\ng4Pv+H11Oh3Tp0/H1taWadOm1Sk/GyXc79b/XcDFixeJjY2lf//+apeiirlz57J8+XIsLMx2mKNG\nkZaWRlFREaGhoQQFBfHBBx9QVFSkdlmqsLe3JywsjK5du9K+fXt+85vftNj/HxWqZ6i3tzc6ne6+\n27Ts/1Eqy83N5bnnnmPlypW0boGTw+7YsYN27drh7+/f4s/ai4qKiI+PZ9y4cURGRhIbG8vXX3+t\ndlmquHbtGqGhoZw9e5akpCSOHTvGzp071S5LVcb8/2iUcA8MDOTcuXOVP8fGxjJgwIDGOLTZ0uv1\njBs3jsmTJzNmzBi1y1HF0aNH2bZtGx4eHkyYMIGDBw8yZcoUtctSRY8ePfDy8uKpp57C3t6eCRMm\nsHv3brXLUoVOp2PAgAH06NGDhx56iPHjx3P48GG1y1JVYGAgcXFxAMTFxREYGHjfbRol3KX/e00G\ng4Hp06fTp0+fFtcToLqlS5eSmppKYmIiGzduZPjw4ayrmA29BfL09CQqKory8nJ27tzJyJEj1S5J\nFcHBwcTExJCdnU1xcTG7d+/mscceU7ssVQUFBREeHk5hYSHh4eF1OjlutGaZiv7vI0eO5KWXXmrR\n/d+PHDnC+vXrOXjwIP7+/vj7+99xdbwlaum9ZT766CNmz55NQEAAdnZ2PN9CJzJo06YNCxYs4Jln\nnmHw4MH4+fkxbNgwtctqNBMmTGDQoEHEx8fTuXNn/v3vfxMaGkpKSgpeXl6kp6czc+bM++5HY2jp\njZ1CCNEMyQVVIYRohiTchRCiGZJwF0KIZkjCXQghmiEJdyGEaIYk3IUQohn6/0A3Dhj0UlGDAAAA\nAElFTkSuQmCC\n"
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In practice, we don't need to implement numerical integration ourselves, as scipy has both basic trapezoid\n",
      "rule integrators and more sophisticated ones. Here we illustrate both:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.integrate import quad, trapz\n",
      "integral, error = quad(f, 1, 9)\n",
      "print \"The integral is:\", integral, \"+/-\", error\n",
      "print \"The trapezoid approximation with\", len(xint), \"points is:\", trapz(yint, xint)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The integral is: 680.0 +/- 7.54951656745e-12\n",
        "The trapezoid approximation with 6 points is: 621.286411141"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "&nbsp;"
    }
   ],
   "metadata": {}
  }
 ]
}