{ "metadata": { "name": "", "signature": "sha256:c277a04292df302bc3cfd5fa01a0f3a93fc22ded39be61d147f33e2fe259db76" }, "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", "\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": [ "Let's first preload the necessary libraries" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "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": 3 }, { "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": 4 }, { "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 = np.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\n" ] } ], "prompt_number": 6 }, { "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": 7 }, { "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": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEmCAYAAAB1S3f/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xu8lGW5//HPxUIERSGPCSJogqCCnITMjIW2CVMjrTR+\nanlC7UBHt8quVyzz11arrf7a/YSdgmmaWJmKGZmnZVqmLgVFQROF5CQQiiKCCFz7j/tZMmtYhzk8\nM888M9/367VezDwz88w1D8O6uO/rPpi7IyIikqtOSQcgIiLposQhIiJ5UeIQEZG8KHGIiEhelDhE\nRCQvShwiIpIXJQ4REcmLEoeIiOSlc9IBiIi0xcz6AqOAg4E/u/vTCYckqMUhIpXtaGAt8DIwIOFY\nJKLEISIVy91/DSwGRgJ3JByORJQ4RKSimFkPM5trZn8xs53dfTFwF9CQcGgSUeIQkUpzCfAg8CHg\n12Z2KPAecEiiUckHTKvjikilMLNOwDJgHLAcGAjsDRwGzHb3FxIMTyJKHCJSMczsKOA37t4n6Vik\nbeqqEpFKMhb4a9JBSPuUOESkkowBnko6CGmfEoeIVAQzqwOOAuYmHYu0TzPHRaRSDAW6A8/GcTIz\n6w8cDgwB7nH3Z+I4r6jFIRXGzJ43s0+08dgvzezycscUp/Y+XzW8X5GOAla5+9qYznciYWTW1cBF\nMZ1TUOKoKmb2jpmtj362mdm7GfcnJh1fLtz9cHf/S1sPRz+p1cHnK4qZLTGzY8vxfq29VwxGA8/H\ndTJ3v8bdnwT6EGafS0zUVVVF3L17820zWwyc6+4PtfZcM+vs7lvKFlx8LNE3r+zr5pTv+pTivT4K\n/CHmcwKcDPyoBOetWWpx1JDof4kXm9lzwHozqzOzS81skZm9bWYvmNlns55/aXT8DTObaWY7Zzze\ny8zuMLPVZvaqmU2Ojp+W0dJZb2bvmdnD0WODzKzRzN6MulFOaiXGY6Pbw8zsmSi2WUDXDj5fQZ8l\nh8+Zfd06tfU5zOwjZrbWzIZlXKM1zd1F0bmOyzr3RWb2XHStZpjZvmY2x8zeMrP7zaxnR5/RzH4F\nHADcE53nouz3a+/aR8/7rpk9a2brzGxW5jXIus47vFdHf68dMbM9gI8A8/N5XQ7n/QzwM6B3nOet\nee6unyr8ITTNj806tgR4hvCPaOfo2OeBD0e3TwXeAfbNeP5z0fM/BDwGXB491gl4Gvg+oeV6IPAK\nMC7rPXcDFgCTgJ2ARcCl0WvGAm8DA7LjBroA/wS+CdQBnwM2Az9s5zMX+lnafKy169bR5wDOA14A\nugH3AT9u6+8luv83wuzoXsCq6L2OiN7rQeAHOX7G1v7Om69nWzH3z/iMfwc+HF2DBcAFuXy/ovO1\n+/eaw/f134BtwJEx/hs4mTC0937ge0n/m6ymn8QD0E+J/mLb/iVyVgevmwuclPH88zMeOx5YFN0e\nDfwz67VTgJkZ9zsRuh7+f3T/GGBl1mt+DUzNjhv4BLA867l/pZ3EUcRnafOx1q5bjp/jbsL/nucB\nO7X19xLdn5hx/3fN1yu6/3Xgzjw+Y1uJo92Yo+f9n4zHrgKm5fL9yuV65PB3dSmwBehWwHf9M8AJ\nwJXA6cCvgIHF/PvRT/s/qnHUnqWZd8zsS8C3gX7Roe7AXm08/zXC/4oB+gK9zOzNjMfrgMxC7I+A\nXYFvRPd7Zb8/oVXRWjdCL8KImOznttmvXsRn6eix7Mdz+Rw3EJLHJHd/v62YI6sybm/Mur+J8DmA\nnD5jW9qKOfNzvp4VR/Y1yPfc+XQPDQP+4e4b83gNZnYAsMDdF5nZDwnJ4y3C36GUiBJH7flgVJKF\n3dV+Qfgf6ePu7mY2l5a/nA/Iur0iur0UWOzurW6uY2ZfBE4jdD1sjQ6vAPqYmXn0X0VCAnqxlVOs\nYMdfPH0JXSKtvV8hn2V5O4+toKXM0Vztfg4z6w5cS0gel5nZ7939TXLXanLM4TO2N+JseXsxt6Kj\n0WuZj+d77tYcQQEzxt39NQAz2xdY7+7rKE2BXTKoOF7bdiX8AvgX0MnMziZMmGpmwFfNrHdUvPwe\nMCt67ElCofhiM+tmodB+uJmNjArD/w2c7C3H5P8deBe42Mx2MrN6wlj7WezocWCLmX0jeu4pwJEx\nf5bbc/icrenoc/w/4El3Px+4F5jezrny0dFnXEUoMLfmiQ5iztbRiKnM93oyz3O3fCOzboStYfPe\nFtbMBprZEcCniVq7ZnZivueR/Chx1DB3XwD8F+GX9OuEX0KPZT6F0Ff9Z0Lh+2Xg/0av3Ur45TAU\neBVYQ/jfcA9Cn3NP4DHbPrLq3qjL5iRCDWEN8HPgTHf/RyuxvQ+cApxF2Dr0VNrZAa6Yz9LBY629\nV5ufw8wmEJYE/0r09O8Awy2/eTSeddtz/IxXAN+PRjZ9J9eY24mhvVbHB+9FqMPkc+5sgwi/iwpZ\no2oc4XtoQFczOxlYXcB5JA8Vu6y6me0KXEfYwKXRwxaSUkbWwVyQNGnvs1TT50wjMzsTmAnsnm+N\nQ5JRyS2OUwjr8p9P+B+siFSnQcCzShrpUdbEYWFi1Sozm591fLyZvWhmL5vZJdHh3mwfqbEVEalW\ng4FHkw5CclfuUVU3EoqmNzcfsLCU8s+BTxJGZzxlZrMJ20f2IUzMquSWUdVy9wOTjiEu7X2Wavqc\nKTUYuCnpICR3Zf2F7O6PAtnDEkcRJlstiQp4s4AJwO+Bz5nZdcDscsYpIuURjWI7AO36lyqVMI8j\ns0sKQktjtLu/C5zT3gvNrDIr+yKSrxVmia5fWXPcveALXgldQEX98k966n01/UydOjXxGKrpR9ez\n458pU6YwefJkXc8y/xSrElocywm1jGZ9CK2OnDQ0NFBfX099fX3ccYlIid15551MmzYNgDfeeIPr\nr7+effbZhyFDhjBixIiEo6s+jY2NNDY2Fn2eSmhxNAH9zayfmXUhLFORc02jOXGISDo8/PDDTJw4\nkRkzZtC9e/cP/v3edNNNjB07ljPPPJOrr7462SCrVH19PQ0NDUWfp9zDcW8jLCE9wMyWmtnZHjbF\n+Tph+ekFwO3uvrCccUmgBBwvXc/WDRgwgKamJq655hpuvfXWD46/+uqr7LfffnTu3Jk33nhjh9fp\nelaOsnZVuXuryy64+xxgTiHnVFdVfHQN46Xr2brevXvz8ssv73B827Zt1NXVAdBaoVzXs3hxdVVV\n7JIjuWi5GKeIpNnPfvYzjjnmGAYNGsRZZ53FrFk5rZEoBTAzvIhRVUocIlIR1q5dy8yZM+nRoweD\nBw/mqKOOSjqkqlXziWPq1KnqqhIRyUFzV9Vll11W24kjzfGLiCSh2BZHJQzHFRGRFFHiEBGRvKQ+\ncTQ0NMQyvExEpNo1NjbGMgFQNQ4RkRqjGoeIiJSVEoeIiOQl9YlDNQ4RkdyoxoFqHCIihVCNQ0RE\nykqJQ0RE8qLEISIieVHiEBGRvKQ+cWhUlYhIbjSqCo2qEhEphEZViYhIWSlxiIhIXpQ4REQkL0oc\nIiKSFyUOERHJS+oTR3vDcVeuhLfeKm88IiKVSsNx6Xg47umnw3HHwTnnlDEoEZEKp+G47Rg1Cp58\nMukoRESqS1UnjtGj4Yknko5CRKS6VHVX1aZNsMce8K9/wS67lDEwEZEKpq6qdnTtCocdBs88k3Qk\nIiLVo6oTB4TuKtU5RETiU/WJY9Qo1TlEROJU9YlDLQ4RkXhVfeLo3x/WrYPVq5OORESkOqQ+cXS0\nkVOnTnDkkWp1iIho5ji5b+T0/e+DGVx+eRmCEhGpcBqOmwPVOURE4lMTLY5Vq2DgQFi7NnRdiYjU\nMrU4crDvvtCjByxalHQkIiLpVxOJAzSfQ0QkLjWTOFTnEBGJR80kDrU4RETiURPFcYANG2DvveHN\nN2HnnUscmIhIBVNxPEe77goDBsC8eUlHIiKSbjWTOEB1DhGRONRU4lCdQ0SkeBWbOMzsQDO7wcx+\nG9c51eIQESlexSYOd1/s7ufFec5Bg+D11+GNN+I8q4hIbSl54jCzmWa2yszmZx0fb2YvmtnLZnZJ\nqeMAqKuDESPgqafK8W4iItWpHC2OG4HxmQfMrA74eXT8UGCimQ0yszPN7Boz61WqYEaPVp1DRKQY\nnUv9Bu7+qJn1yzo8Cljk7ksAzGwWMMHdrwR+FR3bA/hPYKiZXeLuV7V2/sy15evr66mvr283nlGj\nYMaMAj6IiEhKNTY2trtvUb7KMgEwShz3uPvg6P7ngU+5+6To/hnAaHefnOd5c54A2Gz5chg6NOwI\naAVPfxERSa+0TgBMbLp6797QpQssXpxUBCIi6ZZU4lgO9Mm43wdYVsiJOto6tjUalisitShVW8e2\n0lXVGXgJOA5YATwJTHT3hXmeN++uKoArrwybO11zTd4vFRFJvYrvqjKz24C/AQPMbKmZne3uW4Cv\nA/cBC4Db800axVCLQ0SkcOUYVTWxjeNzgDnFnr+hoSGn0VSZRowIix2+/z7stFOxEYiIpENco6tq\nZln1bIcfDjffDMOHxxyUiEiFq/iuqkqlBQ9FRAqT+sRRyKgqUJ1DRGpPqkZVlUoxXVVz58Lpp8OC\nBTEHJSJS4YrtqqrZxLFlC/TsGWaS9+gRc2AiIhVMNY4Cde4Mw4ZBU1PSkYiIpEvqE0ehNQ4IBXLV\nOUSkVqjGQXFdVQC//S3ccAPcd1+MQYmIVDjVOIqIf+NG6NcPGhvD7oAiIrVANY4idOsGX/86/PSn\nSUciIpIeqU8cxdQ4AL76VbjzzjC6SkSkmqnGQfFdVc2+9a2wR8ePfxxDUCIiFU41jhji/+c/w5pV\nr76qOR0iUv1U44hB375w/PEwfXrSkYiIVD61OCLPPQfjx4ctZXfeOZZTiohUJLU4YjJkCBxxBNxy\nS9KRiIhUttQnjmJHVWW6+GL4yU9g27ZYTiciUlE0qop4u6oA3MMyJN/7Hnz2s7GdVkSkoqirKkZm\ncMklcNVVIYmIiMiOlDiynHwyrFkDjz2WdCQiIpVJiSNLXR1cdJEmA4qItEU1jlZs3AgHHggPPgiH\nHRb76UVEEqUaRwl06waTJ4cRViIi0lLnpAMoVkNDA/X19dTX18d63q98BQ4+GJYtg/33j/XUIiKJ\naGxsjGX6grqq2vGd70CnTlp2XUSqixY5LGH8r70GQ4fCK6/Ahz5UsrcRESkr1ThK6IAD4MQTtfih\niEgmtTg6MH8+jBsXFj/s2rWkbyUiUhZqcZTY4MFhr46bb046EhGRyqAWRw4eeQQmTYKFC8MEQRGR\nNFOLoww+8QnYYw+4++6kIxERSZ4SRw7MwpLrWvxQRKQKEkec+3G0Z8IEePNN+MtfSv5WIiIlof04\nKF+No9n118Ndd8G995btLUVEYqcJgGWMf9OmsPjhn/8cRluJiKSRiuNl1LUrfOMbWvxQRGqbWhx5\nWrcODjoI5s0LM8tFRNJGLY4y69kTzjkHrrkm6UhERJKhFkcBli2DIUNg0aIwv0NEJE1K2uIws05m\n9rFCT16t9t8/DM+dNi3pSEREyq/DFoeZzXP3oWWKJy9JtTgAXngBjjsuLH7YrVsiIYiIFKQcNY4H\nzOzzZlbwm1Sjww6DI4+Em25KOhIRkfLKpcXxDrALsBXYFB12d9+9xLF1KMkWB8Bjj8FZZ8FLL2nx\nQxFJj5K3ONy9u7t3cved3H236CfxpFEJjj4a9tkHZs5MOhIRkfLpnMuTzGwC8AnAgUfc/Z6SRtXy\nfU8AdgdmuPv95XjfXJmF3QHHj4dt2+CCC5KOSESk9HLpqroSOBK4FTDgi0CTu08pfXgfxNAT+Km7\nn5d1PNGuqmaLFsGnPgVnnw3f+15IKCIilarka1WZ2XxgqLtvje7XAfPcPefVmsxsJqHlsDrzdWY2\nHrgWqANucPer2nj9T4Fb3H1e1vGKSBwAK1eGlseYMXDttdBJUytFpEKVY1SVAz0z7veMjuXjRmB8\n5oEoAf08On4oMNHMBpnZmWZ2jZn1suAqYE520qg0++0XdgqcNw/OOAM2b046IhGR0sglcVwBPGNm\nvzSzm4Cngf/M503c/VHgzazDo4BF7r7E3d8HZgET3P1X7v5td18BTAaOAz5vZhVfQejZE+67DzZs\ngJNOgnfeSToiEZH4tVscN7NOwDbgKEKdw4FL3X1lDO/dG1iacX8ZMDrzCe7+M+Bn7Z0kc1OS+vp6\n6uvrYwitcN26wR13wPnnhwmC994Le+2VaEgiUuMaGxtj3fAulxrH0+4+oug3MusH3NNc4zCzzwHj\n3X1SdP8MYLS7T87jnBVT48jmDlOmhH3K77tPK+mKSOUoR43jfjO7yMz6mNkezT+FvmGG5UCfjPt9\nCK2OqmAGV14J550HH/84LFyYdEQiIvHIZR7HFwldVF/LOObAQUW+dxPQP2qJrABOAybme5KGhoaK\n6KJqy3e/C3vvDWPHhtbH6NEdv0ZEpBTi6rJqt6sqqnF8wd1vL+pNzG4DxgB7AquBH7j7jWZ2PNuH\n485w9yvyPG/FdlVlu/fesDzJLbeEOR8iIkkpxzyOWGocpZCmxAHw17/CKaeEeR4T825biYjEo9jE\nkUtX1f1mdhFwO7Ch+aC7v1Hom8ap0ruqMh19NDz4IBx/PKxZE/YvFxEpl7J0VQGY2RJamfDn7gcW\n/e5FSluLo9mSJTBuHJx2Gvzwh1qiRETKq+RdVZUsrYkDYPVq+PSnYcQIuO46LcsuIuVTsuG4ZnZx\nxu0vZD2W18zxUmpoaIh1Yku57LMPPPwwvPIKnHoqbNrU8WtERIrR2NjYYtJ0odpscZjZXHcfln27\ntftJSXOLo9l774W1rdauhbvugt2104mIlFg5JgBKCe28M8yaBYccAvX1sGpV0hGJiLRPiaMC1NWF\nOsdnPhNmmS9enHREIiJta2847hAzWx/d7pZxG6BbCWPKS5qG47bHDBoawoKIxxwDf/wjDBmSdFQi\nUk3KNhy3klVDjaM1s2aFOR533BGSiIhInFTjqEJf/GJYmuSUU+CesuzuLiKSOyWOCjVuXFjfatIk\n+OUvk45GRGS7XJYcqWjVUuNozahR0NgYFkVcswb+/d+TjkhE0kw1Dqq3xpFt2bLQAjnhBPjxj7VE\niYgUR0uOpDj+fKxdCyeeGOZ73HADdE59W1FEkqLieI3Yc0944IEwQfDkk+Hdd5OOSERqlRJHiuy6\nK8yeDT16hLrHm28mHZGI1KLUJ460LnJYqJ12gptvDqvqjhkDK1YkHZGIpEXJFzlMg1qqcWRzhyuv\nhOuvh/vug/79k45IRNKiHDsASgUygylTYO+9Q8vjD3+A4cOTjkpEaoESR8qdd14onI8fD7ffDmPH\nJh2RiFS71Nc4JIyy+s1vwla0v/990tGISLVTi6NK1NeHWscJJ8C//gXnn590RCJSrZQ4qsiwYfCX\nv2xfouQ//kOzzEUkfqnvqqq14bgdOfhgeOyx0HX1rW/Btm1JRyQilULDcant4bgdWbcu7Ci4//5h\ndd0uXZKOSEQqhZYckVb17BlqHhs2hASyYUPSEYlItVDiqGLduoVdBHv1guOOCwsliogUS4mjynXu\nDDNmhFFXxxwDS5cmHZGIpJ1GVdUAs7A8yd57w8c/Dn/6EwwalHRUIpJWShw15LvfDclj7Fi4+24Y\nPTrpiEQkjdRVVWO+9KXQdXXiiaF4LiKSLyWOGnTCCXDXXSGJ3HZb0tGISNqoq6pGHX00PPggHH98\nWKJk8uSkIxKRtEh94mhoaKC+vp76+vqkQ0mdww+HRx+FcePCEiWXXaYlSkSqWWNjYywrbWjmuLB6\nNXz609C3L1x0EXz0o0ogItVMM8elaPvsA42N8LGPhbrHsGEwfTqsX590ZCJSidTikBa2bYOHHgqJ\n46GH4NRT4StfgSOOSDoyEYmLWhwSq06d4JOfhN/9Dp5/Hnr3DkN3jzoKbroJNm5MOkIRSZpaHNKh\nLVtgzhyYNg2efDJ0Z11wARxySNKRiUgh1OKQkuvcGU46Cf74R3jqKejaFcaMgWOPDft+bN6cdIQi\nUk5qcUhBNm+GO+8MtZAXX4RzzoFJk6Bfv6QjE5GOqMUhiejSBU47DR5+OPy8+y6MHBnqIX/4A2zd\nmnSEIlIqanFIbN59N3RdTZ8OK1eGFsi558J++yUdmYhkUotDKsYuu8BZZ8Hf/x7Wwlq6FA49FL7w\nhbC8ifY/F6kOanFISb39Ntx6axiRtWlTGI111lmw555JRyZSu4ptcVRs4jCzgcA3gT2B+9x9RivP\nUeJICXd4/PHQjTV7dtgH/cILw/wQLW8iUl5VmziamVknYJa7n9rKY0ocKbR2LfzylyGJ7LJLSCCn\nnw677550ZCK1oeJrHGY208xWmdn8rOPjzexFM3vZzC5p47UnAfcCs0odp5TPnnuG3Qhfegmuvjos\nbdK3b+jGmjs36ehEpCMlb3GY2THAO8DN7j44OlYHvAR8ElgOPAVMBEYCw4GfuPuKjHPc7e4TWjm3\nWhxVYuXKsDPhL34BvXqFVsipp4YWiYjEKxVdVWbWD7gnI3EcBUx19/HR/UsB3P3KjNeMAU4BugIL\n3f3aVs6rxFFltm4NM9SnT4cnnoAzzwwtkYEDk45MpHoUmziS2sipN7A04/4yYHTmE9z9EeCRjk7U\n0NDwwW1t6JR+dXVheZOTToIlS+D666G+HgYNCq2Qk08Okw9FJHdxbeDULKkWx+eA8e4+Kbp/BjDa\n3fPawFQtjtqweXOYFzJ9OixYsH15kwMPTDoykXSq+OJ4G5YDfTLu9yG0OkR20KVLqHc89FDYcGrT\nJjjyyLBr4ezZWt5EpNySShxNQH8z62dmXYDTgNmFnKihoSHWJphUtoEDw0ispUvDWllXXBFaHpdf\nDitWdPx6kVrW2NjYonu/UOUYVXUbMIYwkW818AN3v9HMjgeuBeqAGe5+RQHnVleVMHcu/M//wO23\nw3HHhVrIsceGTalEZEepGFVVKmbmU6dOVVFcgLC8ya9/HZY3effdMBrr7LO1vIlIs+Yi+WWXXVbb\niSPN8UtpuIeFFqdPh7vvDiO0LrwQPvYxLW8iAmpxKHFIu9auDXulT58edi688EI44wwtbyK1TYkj\nxfFL+biHDaemTYMHHghLvV94IQwfnnRkIuWX1uG4sdGoKsmFWSiY//a3YS5I375hMuHo0XDjjaEm\nIlLtUjOqqpTU4pBibN0Kc+aEbqzHHw9dWBdeGGapi1Szmm9xiBSqrm77HunPPAPdu8PYsWGJk1mz\n4L33ko5QpDKlPnGoq0ri0Lcv/OhH8Npr8LWvhTWyDjgApkyBxYuTjk4kHuqqQl1VUlovvRQmFt58\nc1ji5MIL4YQToHNSS4OKxESjqlIcv6TDxo2hqD5tGixbFhZYPPdc6N076chECqMah0iJdesGX/pS\nKKDfc0/YdOrww+GUU+D++2HbtqQjFCkvtThECrB+Pdx6a2iFbNiwfXmTvfZKOjKRjtV8i0PFcUnC\nbruFmse8eXDLLfDCC3DwwWFI72OPhQmHIpVGxXHU4pDK8sYb25c32WknmDABRoyAkSOhTx+tkyWV\nQ8XxFMcv1ckdHnkkLHHy9NPw1FPh2MiR4ac5mfTqpWQiyVDiSHH8UhvcYflyaGoKiaSpKfx07twy\nkYwcCR/+cNLRSi1Q4khx/FK73MNkw8xE0tQURnA1J5HmpLLPPklHK9Wm5hOHNnKSauEOS5a0TCRP\nPx2WgM9OJtqcSgqhjZxQi0Oq37Zt8OqrLRPJM8/AHnu0TCQjRsCHPpR0tJIWNd/iSHP8IoXYtg1e\nfrllzWTuXNh335Y1k+HDoUePpKOVSqTEkeL4ReKydWtYWyuzZvLss2FZlMzi+7BhYQ6K1DYljhTH\nL1JKW7bAiy+2rJnMnx9W/c2smQwdCrvumnS0Uk5KHCmOX6Tc3n8/7ICYWTN5/nk46KCWNZMjjoBd\ndkk6WikVJY4Uxy9SCTZvDskjs2aycCH079+yZjJkCHTtmnS0EodiE0fqdxZoaGjQcFyRInTpEgrp\nw4dvP7ZpU+jWak4kv/gF/OMfMHBgy5rJ4MHh9ZIOzcNxi6UWh4jkZONGeO65ljWTV16BQw9tWTM5\n7LCwVpdULnVVpTh+kbTbsCGM3sqsmSxZEvYryayZHHqodk6sJEocKY5fpBq9806YV5JZM1m6NBTc\nM2smAwdCXV3S0dYmJY4Uxy9SK95+O8x4z5xn8vrrYShwZs1kwADolPpdgiqfEkeK4xepZevWhWSS\n2c21Zk2YpJhZM/nIR5RM4qbEkeL4RaSltWt3TCbr1oURX5k1k4MO0l4mxVDiSHH8ItKxNWu2d3E1\n/7lhw/YuruY/+/ZVMsmVEkeK4xeRwrz+ekgimTWTzZt33Bhr//2VTFqjxJHi+EUkPitW7LgxFrRM\nJM1b9tY6zRzXzHERISSEXr3gpJPC/cwte5ua4LrrQmJp3rI3s2ZSK1v2auY4anGISH6at+zN3v99\nl11aJpKRI2HvvZOOtnTUVZXi+EUkee6weHHLZPL002ETrMxEUk1b9ipxpDh+EalM2Vv2NjWFYcJ7\n7bXjLotp3LJXiSPF8YtIemRu2dvcKpk7N9RHMgvww4fD7rsnHW37lDhSHL+IpFvzlr2Z3VzNW/Zm\nFuCHDYPu3ZOOdjsljhTHLyLVZ8uWsBFWZvF9/vwwQTGzAJ/klr1KHCmOX0RqQ/aWvU1N8MILYR2u\nzJrJEUdAt26lj0eJI8Xxi0jtytyyt7mra+HCsEJwZs1kyBDYeed431uJI8Xxi4hkat6yN7Nm0rxl\nb2bN5PDDi9uyV4kjxfGLiHRk48ZQcM+smbzyStiiN7Nmks+WvUocKY5fRKQQ2Vv2NjXBP/8Jgwe3\nrJkMGtT6lr1VnTjMbFegEWhw93tbeVyJQ0QEWL8e5s1rWTNZtiwU3DNrJoccAp07V3fiuAxYDyxU\n4ii9xsZGLRYZI13PeOl65u+tt3bc//311+Gdd4pLHCXfkNHMZprZKjObn3V8vJm9aGYvm9klrbzu\n34AFwJpSxyhBHKtmyna6nvHS9cxfjx5QXw8XXQS33RZmvr/2WvHnLcey6jcC/w3c3HzAzOqAnwOf\nBJYDT5n1NY1EAAAE+ElEQVTZbGAkMBz4CTAG2BU4FNhoZn9U80JEpDhxrK1V8sTh7o+aWb+sw6OA\nRe6+BMDMZgET3P1K4FfRc74fPfZlYI2ShohIZShLjSNKHPe4++Do/ueBT7n7pOj+GcBod5+c53mV\nTERECpDGHQBj+YVfzAcXEZHClLw43oblQJ+M+32AZQnFIiIieUgqcTQB/c2sn5l1AU4DZicUi4iI\n5KEcw3FvA/4GDDCzpWZ2trtvAb4O3EcYcnu7uy8sdSwiIlK8kicOd5/o7r3cfWd37+PuN0bH57j7\nIe5+sLtfke95O5oHIh0zsyVm9pyZzTWzJ6Nje5jZ/Wb2DzP7s5n1TDrOStTa/KT2rp2ZTYm+qy+a\n2bhkoq5cbVzPBjNbFn0/55rZ8RmP6Xq2wcz6mNnDZvaCmT1vZt+Ijsf2/Uyqq6ooGfNAxhPmeUw0\ns0HJRpVKDtS7+zB3HxUduxS4390HAA9G92VHNxK+f5lavXZmdiihO/bQ6DXXmVkq/+2VUGvX04Gr\no+/nMHefA7qeOXgf+La7HwZ8FPha9Psxtu9nWi/2B/NA3P19YBYwIeGY0ip7ZNpngJui2zcBny1v\nOOng7o8Cb2YdbuvaTQBuc/f3o7lLiwjfYYm0cT1hx+8n6Hq2y91fd/d50e13gIVAb2L8fqY1cfQG\nlmbcXxYdk/w48ICZNZnZpOjYvu6+Krq9Ctg3mdBSqa1r14uWowb1fc3dZDN71sxmZHSt6HrmKJpD\nNwx4ghi/n2lNHJr4F4+j3X0YcDyhOXtM5oPRbH1d6wLkcO10XTs2DTgQGAqsBP6rnefqemYxs+7A\nHcA33X195mPFfj/Tmjg0DyQG7r4y+nMNcCehebrKzD4MYGb7AauTizB12rp22d/X/aNj0g53X+0R\n4Aa2d5/oenbAzHYiJI1fuftd0eHYvp9pTRyaB1IkM9vFzHaLbu8KjAPmE67jl6OnfRm4q/UzSCva\nunazgS+aWRczOxDoDzyZQHypEv1ya3Yy4fsJup7tMjMDZgAL3P3ajIdi+34mteRIUdx9i5k1zwOp\nA2ZoHkje9gXuDN8xOgO3uvufzawJ+I2ZnQssAU5NLsTKFc1PGgPsZWZLgR8AV9LKtXP3BWb2G8Kc\npS3AV7VoZ0utXM+pQL2ZDSV0mywGLgBdzxwcDZwBPGdmc6NjU4jx+1nRGzmJiEjlSWtXlYiIJESJ\nQ0RE8qLEISIieVHiEBGRvChxiIhIXpQ4REQkL0ocIiVmZmPM7Kik4xCJixKHSOmNBT6WdBAicdEE\nQJE8RSuOzgEeJSSE5YSlqf8EzCPMgO4MnAOsAR4Htka3JwP7EWaabwXecvcxZf0AIkVK5ZIjIhXg\nYOA0dz/fzG4HPkdYGqObuw+LVhqe6e6DzWw6sN7drwYws+eAce6+0sx2T+wTiBRIXVUihVns7s9F\nt58G+kW3b4MPNiba3cx6RMczNyT6K3CTmZ2H/vMmKaTEIVKY9zJub6XtBLBDX7C7fwX4PmEp66fN\nbI/4wxMpHSUOkXidBmBmHwfWufvbwHpgt+YnmNlH3P1Jd59KqHvsn0ikIgVSM1mkMG2NKtlkZs+w\nvTgOcA/wOzP7DPAN4Ntm1p/QffVARpeXSCpoVJVITMzsYeC77v5M0rGIlJK6qkREJC9qcYiISF7U\n4hARkbwocYiISF6UOEREJC9KHCIikhclDhERycv/Atvsm8m9Y+A3AAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 9 }, { "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 = np.linspace(0, 10, 200)\n", "y = f(x)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 10 }, { "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": 11 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot both the function and the area below it in the trapezoid approximation" ] }, { "cell_type": "code", "collapsed": false, "input": [ "plt.plot(x, y, lw=2)\n", "plt.axis([0, 10, 0, 140])\n", "plt.fill_between(xint, 0, yint, facecolor='gray', alpha=0.4)\n", "plt.text(0.5 * (a + b), 30,r\"$\\int_a^b f(x)dx$\", horizontalalignment='center', fontsize=20);" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VOXZ+PHvnWVCSEgCBAIJIPsqEJCt2Na4o1JFWxde\ntGjdqliRt1rE/hSoShFcUCu0taK81bIpRVBcEIwbCLIphIQQCCEBEiCBkHWyPb8/zmSDEGAyyZkk\n9+e65pozZ72ZC26euc9znkeMMSillGpafOwOQCmllOdpcldKqSZIk7tSSjVBmtyVUqoJ0uSulFJN\nkCZ3pZRqgmpN7iKyUEQyRGRnDdv+KCJlItKmyrppIrJXRBJE5Jr6CFgppdS5navl/jYw5vSVItIZ\nuBpIqbKuP3A70N91zHwR0V8GSillg1qTrzHmG+BEDZteBv502rqbgMXGmGJjzAEgCRjhiSCVUkpd\nmAtuWYvITUCaMean0zZFAmlVPqcBUXWITSmllJv8LmRnEWkJPIVVkqlYXcshOraBUkrZ4IKSO9AD\n6Ar8KCIAnYCtIjISOAR0rrJvJ9e6akREE75SSrnBGFNbY7qaCyrLGGN2GmMijDHdjDHdsEovQ40x\nGcAq4A4RcYhIN6AXsPks59GXMUyfPt32GLzlpd+FfheN4bs4dMjg52fw8TEkJzfstS/UubpCLgY2\nAL1FJFVE7jk9T1dJ2LuBZcBu4BPgYeNOREop5aXeeANKSuDmm6FrV7ujqV2tZRljzPhzbO9+2udZ\nwCwPxKWUUl4lLw/+/ndrecoUe2M5H9oP3UYxMTF2h+A19LuopN9FJW/6LhYuhKwsGDECRo+2O5pz\nk4aunIiIVmuUUo1KcTH07AkHD8IHH8AttzR8DCKCqa8bqkop1RwtXWol9j59YNw4u6M5P5rclVKq\nFsbAnDnW8hNPgE8jyZpallFKqVqsWQM33ACRkbB/PwQE2BOHlmWUUsqDXnjBen/sMfsSuzu05a6U\nUmexcaPVMyY01Kq5h4TYF4u23JVSykOee856f/hhexO7O7TlrpRSNdiyBYYPh6AgOHAAwsPtjUdb\n7kop5QHPPmu9P/yw/YndHdpyV0qp02zfDkOHQmAgJCdDRITdEWnLXSml6qy81v7733tHYneHttyV\nUqqKnTth0CCr22NyMnTsaHdEFm25K6VUHZTX2h94wHsSuzu05a6UUi47dsCQIVarPSkJOnWyO6JK\n2nJXSik3Pf209f7QQ96V2N2hLXellKLyadSgIGsMmfbt7Y6oOm25K6WUG/78Z+v9sce8L7G7Q1vu\nSqlmb906uOoqCAuzesiEhdkd0Zm05a6UUhfAmMpW+xNPeGdid4e23JVSzdqKFfDrX1ulmH37IDjY\n7ohqpi13pZQ6T0VFMHWqtTx9uvcmdndocldKNVv/+IfVn713b7j/fruj8axak7uILBSRDBHZWWXd\nXBGJF5EfRWSFiIRW2TZNRPaKSIKIXFOfgSulVF1kZ8PMmdbynDng729vPJ52rpb728CY09Z9Dgww\nxgwGEoFpACLSH7gd6O86Zr6I6C8DpZRXmj0bMjPhF7+AG2+0OxrPqzX5GmO+AU6ctm6tMabM9XET\nUP4c103AYmNMsTHmAJAEjPBsuEopVXcHD8Irr1jLL74Ict63KRuPurasfwescS1HAmlVtqUBUXU8\nv1JKedyTT4LTCXfcASOaaBPUz90DReTPQJEx5j+17FZjn8cZM2ZULMfExBATE+NuGEopdUG++QYW\nL4YWLazSjLeKjY0lNjbW7ePP2c9dRLoCq40xA6usuxu4H7jSGFPoWvckgDFmtuvzp8B0Y8ym086n\n/dyVUrYoLYVLLoEff4QZM6zuj41FvfdzF5ExwBPATeWJ3WUVcIeIOESkG9AL2Hyh51dKqfry5ptW\nYu/SxXoatSmrtSwjIouBy4BwEUkFpmP1jnEAa8W6C7HRGPOwMWa3iCwDdgMlwMPaRFdKeYusrMph\nBl56CVq2tDee+qbDDyilmoVJk2D+fLj8cmugsMbWQ+ZCyzKa3JVSTd4PP8DIkeDjA9u3w8CB5z7G\n2+jYMkopVUVJiTUfqjEwZUrjTOzu0OSulGrS5s0rZccO6yZqlV7YTZ6WZZRSTda+fcUMGCA4nX6s\nXg1jx9odkfu0LKOUUoDT6eT22zNwOv0YPDixUSd2d2hyV0o1Ofn5+Tz++Pds3dqJli1LGTcu1u6Q\nGpwmd6VUk5Kbm8ubb65k0aKRADzyyEFCQ/NsjqrhaXJXSjUZ2dnZLF68mKVLf05OTguGDz/FuHFH\n7Q7LFm4PHKaUUt4kKyuLZcuWsW9fNBs3diEwsJSnn07Bp5k2YZvpH1sp1ZQcPXqUxYsX07LlRbz1\n1iUATJ6cRmRkkc2R2UeTu1KqUTty5AhLlizhoou68s47ozhxwp9hw05xyy3H7Q7NVlqWUUo1Wqmp\nqaxYsYI+ffrw3Xd9+fLL1gQFNe9yTDlN7kqpRmn//v2sWrWKAQMGkJPTkRdf7AzA1KkHiYpqvuWY\ncprclVKNzp49e1izZg2DBg0iKKg1jz3WjcJCX8aMyeT667PsDs8raHJXSjUqu3bt4osvvmDo0KEE\nBwfzxhsd2b07iI4dnTz55EG7w/MamtyVUo3Gtm3b+Prrrxk6dChBQUFs2BDCO+90wMfH8OyzyQQH\nl9kdotfQ5K6UahQ2btzI5s2bGT58OC1atCA93Z+nn+6GMcKDDx4iOrr5PYVaG03uSimvZozhm2++\nYceOHQwfPpyAgACKi4Vp07qTne3H6NHZ/O536XaH6XU0uSulvJYxhvXr1xMfH8+IESPw9/cH4LXX\noti5M5iIiCL+8pfkZt/tsSaa3JVSXqmsrIzPPvuMAwcOMHz4cPz8rHT1+eetWbw4Al9fw1//up+w\nsFKbI/VOmtyVUl6npKSENWvWcOTIES655BJ8fX0BiI9vycyZXQGYMiWVQYO0zn42mtyVUl6luLiY\nDz/8kJMnTzJ06FB8XDWX48f9ePzxHjidPtx003Fuv/2YzZF6N03uSimv4XQ6+e9//0thYSHR0dGI\nWLPKFRUJf/pTDzIyHAwenMvUqQeR855wrnmq9TaEiCwUkQwR2VllXRsRWSsiiSLyuYiEVdk2TUT2\nikiCiFxTn4ErpZqWgoICli9fTnFxMQMHDqxI7MbAc89dxE8/WTdQ58zZh8Oh8zCfy7nuMb8NjDlt\n3ZPAWmNMb2Cd6zMi0h+4HejvOma+iOg9bKXUOeXm5rJkyRJ8fHzo379/tW0LFkSyZk1bAgNLeeml\nJNq2LbEpysal1uRrjPkGOHHa6huBRa7lRcA41/JNwGJjTLEx5gCQBIzwXKhKqaYoOzubJUuW0LJl\nS/r06VNt2wcfhLNwYUd8fQ2zZ++nb98Cm6JsfNxpWUcYYzJcyxlAhGs5Ekirsl8aEFWH2JRSTVxW\nVhaLFy+mdevW9OzZs9q2r78O5YUXugDw1FMpXHrpKTtCbLTqdEPVGGNEpLbiV43bZsyYUbEcExND\nTExMXcJQSjVCx44dY/ny5URFRdGpU6dq27ZuDWbatO6UlQkPPHCYm27KtClK+8TGxhIbG+v28e4k\n9wwR6WCMSReRjkD57LOHgM5V9uvkWneGqsldKdX8HDlyhPfff59u3brRsWPHatt27gxiypSeOJ0+\n3HzzMe6//4hNUdrr9IbvzJkzL+h4d8oyq4CJruWJwMoq6+8QEYeIdAN6AZvdOL9SqglLTU1l2bJl\n9OzZ84zEvmdPII8+2pP8fF+uuy6TJ5/ULo/uqrXlLiKLgcuAcBFJBZ4BZgPLRORe4ABwG4AxZreI\nLAN2AyXAw8YY7a+klKqQnJzMhx9+SP/+/Wnbtm21bUlJLZg0qRc5OX5cfvkJpk8/gOvBVOWGWpO7\nMWb8WTZddZb9ZwGz6hqUUqrp2bNnD5988gkDBw4kLCys2raEhEAmTepdMcrjrFnJ+OkjlnWiX59S\nqt7FxcWxdu1aoqOjadWqVbVtP/0UxKOP9iQ314+f//wkL7ywH39//dFfV5rclVL1avv27Xz11VcV\nsydVtWVLMFOm9KSgwJcrrjjB888na2L3EE3uSql68/3337Np0yaGDRtGYGBgtW1r17bmmWe6Ulzs\nw3XXZTJ9+gEtxXiQfpVKKY8zxvDtt9+yffv2itmTKrfBe++1Z948q+f0rbce5fHHU/XmqYdpcldK\nedTZZk8CKCmBV17pzNKl7QGYPDmNO+/M0O6O9UCTu1LKY8pnT0pOTq42exLAyZO+PPVUdzZvDsHf\nv4wZMw5w7bWnD12lPEWTu1LKI0pLS1mzZg2HDx9m2LBhFbMnAezdG8jjj/fg0KEA2rQp5oUX9jNk\nSK6N0TZ9mtyVUnVWXFzMqlWryMrKqjZ7EsCaNW2YNasLhYW+9O+fx5w5++jQodjGaJsHTe5KqTpx\nOp2sXLmS/Px8hgwZUjHJRl6eD3PmdOHjj60nUW+4IZNp01Jo0UK7OjYETe5KKbcVFBSwYsUKSktL\nGTRoUMX6+PiW/PnP3Th4sAUBAWU88cRBbropU2+cNiBN7kopt+Tl5bF8+XL8/f0ZMGAAYM11+uab\nHfm//+tAaanQq1c+s2Yl061boc3RNj+a3JVSF+zUqVMsW7aMVq1aVUyysWtXS2bO7EpyciAihvHj\nM3jkkUMEBGgZxg6a3JVSF+TEiRMsW7aM8PBwunbtysmTvsyfH8V//xuOMUKXLoU888wBoqPz7A61\nWdPkrjyitBTS0iA9HbKyIDPTemVlQX4+FBdbD7CUv4tAYCC0aFH5Cg6GNm2gbVvr1aYNRERASIjd\nfzpVrnz2pMjISCIjO7N8eTsWLIjk1Ck/fH0N//M/6Tz44GG9aeoFNLmrC3LyJOzYYb0SE2H/fti3\nD1JSrMRdH0JCDJ07C507Q+fO0KWL9d69O/TpA+3aoTfqGkB6ejrvv/8+F13Ulfj4fvzv/0aSnGyN\nFzNixCkefzyV7t21tu4tNLmrszIG4uPhm2+s18aNVjI/m44dITKyesu7bVurRe7nB35+ZeTn55CT\nk0lGxnGOH8/Fx6clPj4t8fcPoajIn8xM4cQJITvbn9xcB9nZgZw65U9cHMTF1XzdsDDo3dt69elj\nvV98sfWuA1F5RlpaGh98sILs7NHMn9+P+HhrdMeoKCeTJ6dx+eUn9T9YL6N/9VU1p07BF1/Axx/D\nJ5/AkdOmr2zRAgYOhCFDoF8/6NHDenXtCi1bnnm+rKwsDh06xL59+0hJSSE42I/IyBAuvTSctm3b\nup5iLAZqngDZGKu0s39/MSkphsOHfTh6tAXHjgVy9Ggrjh1ry8mTDjZvhs2nTeoYEAD9+8OgQVbM\ngwZZr4gIT3xTzcfevQeYMSOBjRvvITnZGos9PLyI++47wk03ZeoQvV5KGnomPBHR2fe8TF4erF4N\nS5ZYCb2oqHJbx47wy1/CL34BP/85DBhQe2s4NzeXQ4cOkZyczL59+ygtLSUkJIQ2bdrQrl07HA6H\nx+I2xpCbm0daWhFJST6kpASQlhbEkSMhHDkSTlZWaI3HtWtXmegHDYLBg63/BKoMXKiAw4dhzpzj\nLFrk4ORJ68ZHeHgREyYc5dZbjzaaunpJSQkbNmxgypQpdodSJyKCMea8fx9pcm+mjIHvv4d//AOW\nL7dueoJVux49Gm64wXoNHFh7PdvpdHLkyBFSUlLYt28f2dnZhISEEBYWRrt27c6YnKGhOJ1O0tPz\niY/3JTGxBSkpIaSmhnH4cDiFhWdmcT8/Q9++wuDBVHs11Va+MYaioiIKCwtxOp04nU4KCwvJySni\nyy8drFzZmg0b2lBWZg0j0K1bAXfdlcGYMVk4HI3r368m9waiyd1e+fmwaBEsWAA7d1auHz0a7rgD\nfvMbq7V+NqWlpWRkZJCamsq+fftIT08nKCiI0NBQ2rVrR2hoaMXj596opKSU/ftLiIvzJTExgOTk\nUA4eDOPYsTBq+nfTvn0Z0dHC4MGVib9PH6gyiq2tjDEVybk8QVddzs/PJz8/n4KCgorl8n18fHzw\n8/OjrCyA/fu7sHVrD7Zv70pBgfXryte3jJ//PItbb81ixIgcqgwX06hocm8gmtztkZkJb7wBr78O\nx49b69q3h9/9Du67z6qb18QYQ1ZWFqmpqezfv5+DBw/icDgIDQ0lPDyc1q1bVxv9r7EqKPAhLk6I\ni/Njz54WJCe3IjW1NYWFZ5aR/P0NffuWMmSID0OG+DB4sPULJzzc/etXTdKnJ+jyJF01QRcUFFBQ\nUEBRURG+vr74+fnh6+uLv78/vr6+FYnb4XDgcDjw9/evsuwgLa0V27aF8N13IfzwQwhOZ2Xm7t07\nn2uuyWLs2EzCw0vc/0N5CU3uDUSTe8M6cQLmzoVXX60svQwfDn/8I9x8M9RUAs/JyalWNzfG0KpV\nK9q2bUu7du2qTb7QlJWVweHDDhISAoiL8yMx0Ur6R4+2qnH/0NASevYso08foUePEi66yEnnzvm0\nb5+Hr29htZZ01QSdn59PcXFxteRcnrB9fHzw9/evSM7+/v4EBARU+3yuX0qZmX7s2dOSxMRAdu0K\nYvv2VmRnV79x0rdvHr/4RTbXXHOiyQ0VoMm9gWhybxh5eVZCnzMHsrOtdWPGwNSpcNll1evoTqeT\nw4cPk5KSQlJSEjk5ObRq1YrWrVvTvn17WtbUDaYZy831ISkpkMTEQOLj/UlMDCQlJZjCwrP/pxcU\nVEjbtvmEhxfSvr2TiIgiOnQoJjzcuLqNGsLCSgkKKrugLoVFRUJ2th/Z2b6cPOlHerqDQ4cCSEsL\nqHjPyjozrnbtioiOzmXUqFOMHn2Kdu2a7hC8zTW5u90VUkSmAXcCZcBO4B4gCFgKXAQcAG4zxpx0\n9xrqwhkDS5fC44/DoUPWuiuvhOefh5Ejrc+lpaWkp6dX1M0zMjIIDg4mNDSU7t27Expacy8TZQkO\nLiM6Oq/a4/XGWC3klJQWpKS04ODBgIrl9HQHeXktyMtrwcGDtZ/b19fQokUZAQFlOBxlOBwGh6MM\nY4SSEqG0FEpLreWcHF8KC89dEgsKKqV373x69y6gX788oqNziYoq0n7pTZxbLXcR6QqsB/oZY5wi\nshRYAwwAjhtj5ojIVKC1MebJ047Vlns92bkT/vAH+Oor6/PQoVbL/YorDJmZmdXq5i1atCAkJITw\n8HDatGlTbXIF5VnGwIkTfmRkOEhPd5CR4V+xbLW6rZZ3drYfBQUXdv/Cz6+M0NBSQkNLCAkpISKi\nmKgoJ506OSve27UrbrQ3Qz1BW+4X5hTWkyctRaQUaAkcBqYBl7n2WQTEAk/WdALlOU4nPPsszJ5t\njfESHg7PPFPA5Zcf4ODBZN54Y39F3Tw8PJzRo0dXm9tS1S8RaNOmhDZtSujXL7/WfYuKBKfTh6Ii\noajIB6fTevfxMfj6Wl02fX0Nfn6GoKBSWra8sDKOaj7c+hdujMkSkZeAg0AB8JkxZq2IRBhjMly7\nZQBNtJew99i82erxEhcHIoabbjpETMw6Skqy2Lw5lNatWxMdHU1gYKDdoarzYJVhSu0OQzUBbiV3\nEekBPAZ0BbKB5SJyZ9V9jDFGRGqsv8yYMaNiOSYmhpiYGHfCaNZKSuDZZw3PPQdlZUJ4+Anuuedr\nLrmkkPbtuxAScrHdISql6iA2NpbY2Fi3j3e35n47cLUx5j7X57uAUcAVwOXGmHQR6Qh8aYzpe9qx\nWnOvo9RUmDDBGsxLxHDLLSlMnnycli3197lSp2uuNXd3b7MkAKNEJFCsTrZXAbuB1cBE1z4TgZVu\nnl+dxapV1lOS33wDISG5vPzyTqZNy9TErpSqxt2a+48i8n/AFqyukNuAfwKtgGUici+urpAeirPZ\nKyuDv/wFZs60Pg8ceJCpUxPo27etvYEppbyS210mjDFzgDmnrc7CasUrDzp1Cu66y2q1+/jAXXfF\nc/nlP9C37wC7Q1NKeSntD+flkpOt0Rnj461JKWbN2k9JyRf07z/K7tCUUl6sGT/a4P22bIFRo6zE\nPmAAfPppJk7nKqKjo/WhI6VUrTRDeKk1a6wxYI4etYYPWL/eyc6d/6Vnz5461otS6pw0uXuhd96B\nG2+0RnH87W+tRL9581ocDgcdaxtsXSmlXDS5e5n58+Gee6xhBP78ZyvRJyT8xIEDB+jXr5/d4Sml\nGglN7l7k5Zdh0iRr+cUX4bnn4PjxY6xfv55BgwZpnV0pdd60t4yXmDXLaqmDNWPSww9b46x/+OGH\n9OjRw7a5SJVSjZMmdy/w0ktWYheBt96yyjIAa9euxc/Pj8jISHsDVEo1Ovo732YLFlgTawAsXFiZ\n2Hfu3ElycjL9+/e3LzilVKOlyd1GixZZ5RewbqTefbe1fOyYVWcfPHhwk5h8WinV8DS52+TDD61x\n2MG6efrQQ9ZyUVERq1atolu3blpnV0q5TZO7DTZuhDvusAYDmz4d/vjHym1ffPEFvr6+REVF2Reg\nUqrR0+TewBIT4Ve/gsJCuO8+K7mX27VrF0lJSVpnV0rVmSb3BpSRAWPGQGYmXH+9dTO1fP7L48eP\ns27dOq2zK6U8QpN7A3E6Ydw4a5THYcNg6VIon6O6vM7etWtXgoOD7Q1UKdUkaHJvAMbA738P338P\nnTvDRx9B1Ry+bt06fHx86NSpk31BKqWaFE3uDWDePGuMmJYtrQk3IiIqt8XFxWmdXSnlcZrc69ln\nn1U+pPTOOxAdXbktMzOTL774gkGDBmmdXSnlUZrc61FycmWXx2eegVtvrdxWXFysdXalVL3R5F5P\nCgvhN7+BkyetsdmrdnkEq84uIlpnV0rVC03u9WTKFNi2Dbp1s4YZqDpa7+7du0lMTNQ6u1Kq3mhy\nrwfvvQd//zsEBMD771sTW5fLzMxk7dq12p9dKVWvNLl7WEICPPCAtfzqqzB0aOU2rbMrpRqK28ld\nRMJE5H0RiReR3SIyUkTaiMhaEUkUkc9FJOzcZ2o6nE4YP96a+3TChMokX279+vUAWmdXStW7urTc\nXwXWGGP6AYOABOBJYK0xpjewzvW52fh//w927IDu3a0hfMuHFgCIj48nISGBAQMG2BegUqrZcCu5\ni0go8AtjzEIAY0yJMSYbuBFY5NptETDOI1E2Al98YQ3d6+tr1dxDQiq3ZWVl8fnnnxMdHa11dqVU\ng3C35d4NOCYib4vINhF5U0SCgAhjTIZrnwwg4uynaDoyM2HiRGt5+nQYNapyW3mdvUuXLlpnV0o1\nGHfnUPUDhgKPGGN+EJF5nFaCMcYYETE1HTxjxoyK5ZiYGGJiYtwMwzs89BAcPgyXXgrTplXfFhsb\nS1lZGV26dLEnOKVUoxQbG0tsbKzbx4sxNebf2g8S6QBsNMZ0c33+OTAN6A5cboxJF5GOwJfGmL6n\nHWvcuaa3Wr4cbrvNGgjsp5+sfu3lEhIS+Pzzzxk5ciR+fjoXuVJ2KCkpYcOGDUyZMsXuUOpERDDG\nyLn3tLhVljHGpAOpItLbteoqIA5YDbgKFEwEVrpz/sbi6NHKOVBffLF6Ys/KyuKzzz5j0KBBmtiV\nUg2uLlnnD8B7IuIA9gH3AL7AMhG5FzgA3FbnCL3YI4/A8eNw5ZXVuz1WrbO3atXKvgCVUs2W28nd\nGPMjMLyGTVe5H07jsXy59QoOhn/9q3q3x6+++orS0lKtsyulbKNPqLohMxMmTbKW586Frl0rtyUk\nJBAXF8fFF19sS2xKKQWa3N3ypz/BsWNw+eXw4IOV60+cOKF1dqWUV9DkfoG++goWLgSHwxocrLwc\nU1JSwurVq+ncuTMhVZ9gUkopG2hyvwBOpzUXKsBTT0Hv3pXbvvrqK4qLi7nooovsCU4pparQ5H4B\n5s61Rn3s3RuerPLI1p49e9i1a5fW2ZVSXkOT+3nauxeee85aLh+rHbTOrpTyTprcz9Njj1llmd/+\n1rqRClad/aOPPiIqKkrr7Eopr6LJ/Tx8/DGsWWON9Dh3buX6r7/+GqfTSdeqfSGVUsoLaHI/B6fT\narUDzJwJ7dtby4mJiezcuZOBAwfaF5xSSp2FJvdzmDcPkpKgX7/KB5dOnjzJp59+ysCBA7XOrpTy\nSprca3H4MDz7rLX86qvg7w+lpaWsXr2aqKgoQkND7Q2wmVm2bBmXXXYZu3btsjsUpbyeJvdaTJ0K\neXlw881w9dXWum+++YbCwkKts9vghhtuICAgQKcqVOo8aHI/iy1b4N13rSdRX3rJWpeUlMSPP/7I\noEGD7A2umdqyZQtDhgxB5LyHtFaq2dLkXgNj4IknrOXJk61x2rOzs1mzZg0XX3yx1tlt8v333yMi\nfPrpp8yaNYukpCS7Q1LKa2lyr8FHH0FsLLRpYw0zUFpaWtGfPSwszO7wmoUlS5Zw5ZVXcuedd5KS\nkgLA5s2bmTBhAmPGjOGXv/wl8+fPtzlKpbyXJvfTlJRYoz4CPP00hIVZdfb8/HytszeQLVu28Mor\nrzBv3jzy8vJ49tlnSU9PxxhT0fX0+PHjnDx50uZIlfJemtxP89Zb1vgx3btbU+jt27dP6+wN7LXX\nXuNnP/sZvXv3xhhDREQE8fHxREdHV+yzadMmRo8ebWOUSnk3Te5V5ObCM89Yy7NnQ0FBNh9//LHW\n2RvQrl27iI+P5+qrryYgIICVK1fy/PPPExQUVDFl4cGDB0lKSuLOO++0OVqlvJcm9yrmzbMmvR45\nEm6+2aqzR0ZGap29AX388ccAZ7TKhw8fjo+PDx999BGLFy9mwYIFtGjRwo4QlWoUxBjTsBcUMQ19\nzfNx4kR5rxhYtw58fGLZu3cvQ4cOtTu0ZuWGG24gODiYpUuX2h2KaiJKSkrYsGEDU6ZMsTuUOhER\njDHn3Q9YW+4uc+daif3KK6Fr1/3s2LFDx41pYAcPHuTo0aPVautKKfdocgcyMqzhBQCeeiqPNWvW\nMHDgQPz9/e0NrJn54YcfAHTSE6U8QJM7MGsW5OfD2LGGI0dWEhERoXV2G2zZsgWAfv362RyJUo1f\nnZK7iPiKyHYRWe363EZE1opIooh8LiJenyEPHrRmVgK45ZZt5Ofn0717d3uDaqa2bduGw+GgW7du\ndoeiVKPoXu9dAAARtklEQVRX15b7ZGA3UH6H9ElgrTGmN7DO9dmrzZoFRUUwdmwu2dnfaJ3dJikp\nKWRlZdGjRw98fX3tDkepRs/t5C4inYDrgX8B5XdwbwQWuZYXAePqFF09S02FhQtBxDB48Eouvvhi\nrbPbZPv27QD06dPH5kiUahrq0nJ/BXgCKKuyLsIYk+FazgAi6nD+evfCC1BcDD/7WQqDBjlo3bq1\n3SE1W9u2bQOgZ8+eNkeiVNPg1mOXIjIWOGqM2S4iMTXtY4wxIlJjh/YZM2ZULMfExBATU+Mp6tWh\nQ/Dmm9by1Vf/QI8ePRo8BlVp586dgHck99LSUrdLQyUlJfo0s/KI2NhYYmNj3T7e3b+Fo4EbReR6\noAUQIiL/BjJEpIMxJl1EOgJHazq4anK3y9y5Vq198OC9jBnT2e5wmrUTJ06QlpaGiNj+n+z69evJ\ny8vjV7/6lVvHv/3224wYMYLBgwd7ODLV3Jze8J05c+YFHe9WWcYY85QxprMxphtwB7DeGHMXsAqY\n6NptIrDSnfPXt/R0+Mc/rB8VDz2UqXV2m/30008AtG7dukG6oKampvLYY4/x2muv8de//pXyJ6a3\nbt3K9u3b3U7sAPfccw8LFy4kOTn5vI955ZVXuOGGGxg+fDhbt251+9pKVeWpfu7l5ZfZwNUikghc\n4frsdebONRQWCsOHpzFsmCZ2u5Un94YoyRQXF/PII49w5ZVXkpmZyYcffkheXh65ubm89tprPPLI\nI3U6v5+fH9OmTWP69OmUlJSc1zFTpkxh4sSJOBwO7a2lPKbOxUFjzFfAV67lLOCqup6zPmVmwoIF\nBhAmTz5ldziKynp7r1696v1aGzdu5PDhwwwdOpTu3btXjGXz+uuvc9111xEQEFDna3To0IEePXrw\n0UcfMW7c+XUY2759O/3798fhcNT5+kpBM3xC9Y03oKDAh4ED0+jbt8DucJq90tJS4uPjgYZpuW/d\nupXWrVsTFRXFgAEDGDFiBAUFBaxcuZLrr7/eY9e5/fbbWbRo0bl3dNmxY4cOUqc8qlkl97w8eO01\na/n663fZG4wC4MCBAxQWFiIi9O7du96vFxcXR//+/aut+/bbb4mMjCQkJMRj1+nduzfZ2dkkJCSc\nc9+0tDSOHz+uyV15VLPqs7VwoVWWGTiwgD590oF2dofU7O3evRsAX1/feh32YcaMGWRlZfHjjz/S\ntWtXHn30UaKiopg6dSqbNm2qdaat+Ph41qxZg4+PD4cPH+bpp59mxYoV5OTkcPToUR588EE6depU\n7RgfHx8GDx7M999/T9++fatt++GHH1ixYgWRkZHk5ORUPJV7eg8bd66rVLlmk9yLi+Gll6zl++7L\nRM57VGRVn8qTe7du3eq1f/iMGTM4dOgQ48aNY9KkSdW6mCUmJnLzzTfXeFxaWhqrVq1i6tSpFee5\n++67mTlzJmVlZdx///307duXCRMmnHHsRRddRGJiYrV1K1euZMGCBbz77ru0a9eO9PR0fv3rX9O/\nf/9qk4/U5bpKQTMqyyxbBikp0Ls3XHFFjt3hKJfy5H5667Y+7NmzB+CM8s/hw4crpvA73Xvvvcej\njz5a8bmgoIDQ0FAGDhxIhw4dmDBhwlm7TrZq1YrDhw9XfE5MTGT27Nn88Y9/pF0761djhw4dCAwM\n5JJLLvHYdZWCZpLcjYE5c6zlP/0JdFwq71BaWkpSUhLQMMP8JiYmEhwcTGRkZLX1ubm5BAcH13jM\nXXfdRWBgYMXnnTt3MmLECAAiIiKYPHnyWWv1oaGh5ObmVnyeP38+QUFBXHnllRXr9u/fT3Z29hn1\n9rpcVyloJsn988/hp58gMhJ0TmXvceDAAYqKihCRBknue/bsqfGmrWv6shqPqfofwYEDBzh27BjD\nhg07r+sZYyrOm5OTw8aNGxk5cmS1oQ22bt2Kj4/PGbNP1eW6SkEzSe4vv2y9/+EP4IFuzMpDyuvR\nfn5+DTIa5N69e2tM7q1atSI7O/ucx2/ZsgV/f/9qN1/T0tLOun92dnZFuSc1NZWysrIzbtxu2bKF\nfv36ERgYyKFDhzxyXaWgGST3uDir5d6yJTzwgN3RqKrKk3vPnj3rfQiI7OxsMjIyakzukZGRNSb3\nwsJCXn311YrS0aZNm+jVq1fFg05lZWX8+9//Pus1T506RVRUFABBQUGAVWOvev5t27ZVlGSWLFni\nkesqBc2gt8y8edb7xInQpo29sajq9u3bB8CAAQPq/VrlN1Nrego2Ojq6xrFgvvvuO95991369euH\nn58fqamp1W68Lly4sNabmsnJyYwcORKALl260KtXr4rWeUlJCS+88ALFxcV06tSJrKws2rj+gtb1\nukpBE0/ux45BeQNn8mR7Y1FnKm+ZNkRyT0hIoFWrVjW23H/2s5/xUnk/2SouueQSxo4dS3x8PHv2\n7OGdd95h9uzZzJo1C39/fy677LKzTuZdUlLCTz/9VNHjRUSYPXs2L7/8MhkZGZSVlXHvvfdyySWX\n8NFHHxEfH1+xb12uq1S5Jp3c//53cDrhhhtAJ/jxLjk5ORw7dgwRaZBElZCQwPDhw/HxObMSOWTI\nEDIzMzl27FhFF0WAsLAwpk+fXm3f8x2uOi4ujoiIiGq/FLp06cK88p+SLp06dWLs2LHV1tXlukqV\na7I1d6fTGkcGYMoUe2NRZyovyYSEhNC1a9d6ucY777zDpEmTAKs/fdUuiFU5HA5uu+02Fi9e7LFr\n/+c//+FO7ZqlbNRkk/vixZCRAYMGwRVX2B2NOt3+/fsBzugC6EmffPIJDoeDvXv34u/vf9bkDjBx\n4kQ2bNjAqVN1Hyn0wIEDpKena11c2apJJndj4PXXreXJk9GhBrxQeXIfMmRIvV3jrrvuIjw8nIUL\nFzJ37txap85r0aIFTz/9NM8999xZ+7yfD6fTydy5c3n++ecR/YunbNQka+6bNsG2bdC2LYwfb3c0\nqibl3SDrs+U+duzYM+rZtRkwYAC33HILS5cu5Y477nDrmm+//TaTJk3SAb2U7Zpkcv/b36z3e++F\nKk9wKy+SlJREYGBgg4wpcyFGjRrFqFGj3D7+97//vQejUcp9Ta4sc/QoLF9ulWIeesjuaFRNjhw5\nQk5ODhdffHGtpRKllPuaXHL/17+gqAjGjoV66oSh6qh85qXTR0JUSnlOk0ruJSWwYIG17OoBp7xQ\nXFwcQMUoh0opz2tSyX31akhLg1694Oqr7Y5Gnc3OnTsJCgpqkCdTlWqumlRyL39o6eGHoYYHEZUX\nKCwsJC4ujpEjR9b4tKhSyjOazL+uxERYt84a/fHuu+2ORp3Nli1bKCoq4rLLLrM7FKWaNLeSu4h0\nFpEvRSRORHaJyKOu9W1EZK2IJIrI5yIS5tlwz+5f/7Le77gDwhrsqupcXnzxRcaPH09JSQkAn376\nKSEhIbU+LaqUqjt3W+7FwBRjzABgFDBJRPoBTwJrjTG9gXWuz/XO6YS337aWdcx277J582YKCwsp\nKysjPT2d9evXM378+IqxyZVS9cOth5iMMelAums5V0TigSjgRqD89/YiIJYGSPAffgjHj1vjyGgH\nDO8yePBg2rRpw6lTp/jLX/5Cly5duFvrZkrVuzrX3EWkKzAE2AREGGMyXJsygIi6nv98/OMf1vsD\nD+g4Mt5m0qRJxMXFMW7cOBwOB6+//jp+fjW3KUpKSliwYAHvv/8+S5YsYcqUKTqdnFJuqtPwAyIS\nDHwATDbG5FQdKMkYY0SkxhGYqo5NHRMTQ0xMjNsx7N0L69dbwwxMmOD2aVQ9CQsL42/l40Gcw1//\n+ld69erFb37zG06ePMk///lPHaNFNVuxsbHExsa6fbzbyV1E/LES+7+NMStdqzNEpIMxJl1EOgJH\nazrWkxMPlN9Ivf12vZHamO3du5e1a9cydepUwBp7pnxuUaWao9MbvjNnzryg493tLSPAW8BuY0zV\nqWVWARNdyxOBlacf60lFRXojtanYtGkT0dHROBwOwLoRO2zYMHJycmyOTKnGyd2a+6XAncDlIrLd\n9RoDzAauFpFE4ArX53qzapU1T+rFF0MdBvJTXiA0NJS2bdsCkJ+fz5dffsngwYP55JNPbI5MqcbJ\n3d4y33L2/xiucj+cC7NwofV+3316I7Wxu/baa9mxYwefffYZRUVFXHvttXz33XdeNySwUo1Fox3P\n/dAh+Owz8PfXG6lNgcPh4Omnn7Y7DKWajEY7/MC//w1lZXDjjRAebnc0SinlXRplcjem8kbqPffY\nG4tSSnmjRpncN260Bgrr2BGuvdbuaJRSyvs0yuRefiP1rrvgLA87KqVUs9boknteHixdai1rSUYp\npWrW6JL7Bx9Abq7Vr117ySmlVM0aXXJftMh611a7UkqdXaNK7mlp8OWXEBBgjSWjlFKqZo0quS9e\nbHWD/NWvIDTU7miUUsp7Nark/t571rs+kaqUUrVrNMl91y748Udo3Rquu87uaJRSyrs1muRe3mq/\n9Var5q6UUursGkVyLyurTO533mlvLEop1Rg0iuT+7beQmgpdusCll9odjVJKeb9Gkdzffdd6nzAB\nfBpFxEopZS+vT5VOJyxfbi1rLxmllDo/Xp/c16yBkychOhoGDLA7GqWUahy8Prlr33allLpwXp3c\nT56E1aut+VHHj7c7GqWUajy8Orl/8AEUFcHll0NUlN3RKKVU4+HVyV37tiullHs8ntxFZIyIJIjI\nXhGZ6u55MjLgq6/A4YCbb/ZkhEop1fR5NLmLiC/wN2AM0B8YLyL93DnXihXWk6nXXANhYZ6M0nts\n2bLF7hC8hn4XlfS7qKTfhfs83XIfASQZYw4YY4qBJcBN7pyovG/7bbd5LDavs3XrVrtD8Br6XVTS\n76KSfhfu83RyjwJSq3xOc627IFVLMjfe6LHYlFKq2fDz8PmMJ05StSRTH5NyiAhZWVls27bN8ye/\nAEeOHLE9Bm+h30Ul/S4qeeK7KCsrw8/P06nO+4kxHsnH1slERgEzjDFjXJ+nAWXGmBeq7OO5Cyql\nVDNijJHz3dfTyd0P2ANcCRwGNgPjjTHxHruIUkqpc/LobxVjTImIPAJ8BvgCb2liV0qphufRlrtS\nSinv0KBPqHrqAafGTkQ6i8iXIhInIrtE5FG7Y7KbiPiKyHYRWW13LHYSkTAReV9E4kVkt+s+VrMk\nItNc/0Z2ish/RKTZTLApIgtFJENEdlZZ10ZE1opIooh8LiK1PgHUYMndkw84NQHFwBRjzABgFDCp\nGX8X5SYDu/FQj6tG7FVgjTGmHzAIaJZlTRHpCtwPDDXGDMQq895hZ0wN7G2sXFnVk8BaY0xvYJ3r\n81k1ZMvdYw84NXbGmHRjzA7Xci7WP+BIe6Oyj4h0Aq4H/gWcd2+ApkZEQoFfGGMWgnUPyxiTbXNY\ndjmF1Qhq6eqo0RI4ZG9IDccY8w1w4rTVNwKLXMuLgHG1naMhk7tHHnBqalwtlCHAJnsjsdUrwBNA\nmd2B2KwbcExE3haRbSLypoi0tDsoOxhjsoCXgINYPe9OGmO+sDcq20UYYzJcyxlARG07N2Ryb+4/\nt88gIsHA+8BkVwu+2RGRscBRY8x2mnGr3cUPGArMN8YMBfI4x0/vpkpEegCPAV2xftUGi4hO2eNi\nrJ4wtebUhkzuh4DOVT53xmq9N0si4g98ALxrjFlpdzw2Gg3cKCLJwGLgChH5P5tjsksakGaM+cH1\n+X2sZN8cDQM2GGMyjTElwAqsvyvNWYaIdAAQkY7A0dp2bsjkvgXoJSJdRcQB3A6sasDrew0REeAt\nYLcxZp7d8djJGPOUMaazMaYb1g2z9caY39odlx2MMelAqoj0dq26CoizMSQ7JQCjRCTQ9e/lKqwb\n7s3ZKmCia3kiUGujsMEGXNAHnKq5FLgT+ElEtrvWTTPGfGpjTN6iuZfv/gC852oA7QPusTkeWxhj\nfnT9gtuCdS9mG/BPe6NqOCKyGLgMCBeRVOAZYDawTETuBQ4AtY6Zqw8xKaVUE+TV0+wppZRyjyZ3\npZRqgjS5K6VUE6TJXSmlmiBN7kop1QRpcldKqSZIk7tSSjVBmtyVUqoJ+v9BNdB9o+xqCwAAAABJ\nRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 12 }, { "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.549516567451064e-12\n", "The trapezoid approximation with 6 points is: 621.286411141\n" ] } ], "prompt_number": 14 } ], "metadata": {} } ] }