{ "metadata": { "name": "", "signature": "sha256:cbf49a9ceac0258773ae0a652e9ac7c13e0d042debbf0115cdc081862b5f7308" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Basic Numerical Integration: the Trapezoid Rule" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A simple illustration of the trapezoid rule for definite integration:\n", "\n", "$$\n", "\\int_{a}^{b} f(x)\\, dx \\approx \\frac{1}{2} \\sum_{k=1}^{N} \\left( x_{k} - x_{k-1} \\right) \\left( f(x_{k}) + f(x_{k-1}) \\right).\n", "$$\n", "
\n", "First, we define a simple function and sample it between 0 and 10 at 200 points" ] }, { "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": 1 }, { "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": 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[np.logical_and(x>=a, x<=b)][::30]\n", "yint = y[np.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": [ "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+/AAAIABJREFUeJzt3XlYlXX+//HnAQ6IIIuiBwEVFZFFNnMpp75hDlqmZmWW\nY+Zoy0zapE41WvMrtUUxa1KbdObbZDHV160aNSPKJRzLLbdURFE22VUQZD1s9++PO7ZEksOB+wDv\nx3Wd69yc7X53rnxx874/9+ejUxRFQQghRIdipXUBQgghzE/CXQghOiAJdyGE6IAk3IUQogOScBdC\niA5Iwl0IITqgJsN99uzZGAwGgoKCrnvu7bffxsrKiry8vNrHli9fzqBBg/Dz8+Pbb781f7VCCCFu\nSpPhPmvWLGJiYq57PC0tjZ07d9KvX7/ax86cOcOmTZs4c+YMMTExzJkzh+rqavNXLIQQ4lc1Ge53\n3HEHrq6u1z3+5z//mTfffLPBY9u2bWPatGno9Xq8vb3x8fHh8OHD5q1WCCHETWl2z33btm14eXkR\nHBzc4PHMzEy8vLxqf/by8iIjI6PlFQohhGg2m+a8uKSkhGXLlrFz587ax5qavUCn05lemRBCCJM1\nK9wTExNJSUkhJCQEgPT0dG655RYOHTqEp6cnaWlpta9NT0/H09Pzus/w8fEhMTGxhWULIUTnMnDg\nQC5cuHDzb1B+RXJysjJkyJBGn/P29lZyc3MVRVGUuLg4JSQkRDEajUpSUpIyYMAApbq6+rr33MQu\nO43FixdrXYLFkO+ijnwXdSztu8jIUBQbG0WxslKU5OS23Xdzs7PJnvu0adMYNWoUCQkJ9OnThw8/\n/LDB8/XbLgEBAUydOpWAgADuuece1q5dK20ZIUSH8t57UFkJ998P3t5aV9O0JtsyGzZsaPLNSUlJ\nDX5+6aWXeOmll1pelRBCWJjiYvjHP9TtBQu0reVmyBWqGgoPD9e6BIsh30Ud+S7qWNJ3sX495OXB\niBEwapTW1fw63c+9nLbboU7X5AgbIYSwNBUV4OMDFy/C55/DAw+0fQ3NzU45chdCiF+xaZMa7IMH\nw+TJWldzcyTchRCiCYoCNRfkv/ACWLWT1JS2jBBCNCE6Gu69Fzw8ICkJ7Oy0qUPaMkIIYUYrVqj3\n8+drF+ymkCN3IYS4gQMH1JExzs5qz93JSbta5MhdCCHM5PXX1fs5c7QNdlPIkbsQQjTiyBEYPhwc\nHCAlBdzctK1HjtyFEMIMXntNvZ8zR/tgN4UcuQshxC8cPw5Dh4K9PSQng8GgdUVy5C6EEC1W02v/\n4x8tI9hNIUfuQghRz6lTEBysDntMTobevbWuSCVH7kII0QI1vfannrKcYDeFHLkLIcTPTpyAsDD1\nqP3CBai3LLTm5MhdCCFM9PLL6v3TT1tWsJtCjtyFEIK6q1EdHNQ5ZHr10rqihuTIXQghTPDXv6r3\n8+dbXrCbQo7chRCd3u7d8NvfgouLOkLGxUXriq4nR+5CCNEMilJ31P7CC5YZ7KaQI3chRKf2xRfw\n4INqKyYxERwdta6ocXLkLoQQN6m8HBYuVLcXL7bcYDeFhLsQotP65z/V8ey+vvDkk1pXY15Nhvvs\n2bMxGAwEBQXVPvbCCy/g7+9PSEgIDzzwAAUFBbXPLV++nEGDBuHn58e3337belULIUQLFRTA0qXq\n9ptvgl6vbT3m1mS4z5o1i5iYmAaPjR07lri4OH766Sd8fX1Zvnw5AGfOnGHTpk2cOXOGmJgY5syZ\nQ3V1detVLoQQLRAZCbm5cMcdMGmS1tWYX5Phfscdd+Dq6trgsYiICKx+Xv575MiRpKenA7Bt2zam\nTZuGXq/H29sbHx8fDh8+3EplCyGE6S5ehHfeUbffegt0Om3raQ0t6rmvX7+e8ePHA5CZmYlXvet1\nvby8yMjIaFl1QgjRChYtAqMRHnkERozQuprWYWPqG9944w1sbW353e9+d8PX6G7w63DJkiW12+Hh\n4YSHh5tahhBCNMu+fbBhA3TporZmLFVsbCyxsbEmv9+kcP/oo4+Ijo5m9+7dtY95enqSlpZW+3N6\nejqenp6Nvr9+uAshRFupqoI//UndXrQI+vXTtp6m/PLAd2nN2d+b1Oy2TExMDCtXrmTbtm106dKl\n9vFJkyaxceNGysvLSU5O5vz584zoqH/vCCHapfffh59+gr591atRO7Imj9ynTZvG3r17uXLlCn36\n9GHp0qUsX76c8vJyIiIiALjttttYu3YtAQEBTJ06lYCAAGxsbFi7du0N2zJCCNHW8vLqphl4+23o\n2lXbelqbTD8ghOgU5s6FtWth9Gh1orD2duzZ3OyUcBdCdHg//ggjR4KVFRw/DvWuy2w3ZG4ZIYSo\np7JSXQ9VUWDBgvYZ7KaQcBdCdGirVlVx4oR6ErUzDdSTtowQosNKTKwgMFCH0WjDl1/ChAlaV2Q6\nacsIIQRgNBp5+OEcjEYbQkIS2nWwm0LCXQjR4ZSUlPD88wc5etSLrl2rmDw5VuuS2pyEuxCiQykq\nKuL997cSFTUSgGeeuYizc7HGVbU9CXchRIdRUFDAhg0b2LTpdgoLuzB8+DUmT76kdVmaMHniMCGE\nsCR5eXls3ryZxMRQDhzoi719FS+/nIpVJz2E7aT/2UKIjuTSpUts2LCBrl378cEHtwAwb146Hh7l\nGlemHQl3IUS7lpWVxcaNG+nXz5uPPrqVq1f1DBt2jQceuKJ1aZqStowQot1KS0vjiy++YPDgwfzw\ngx/ffeeKg0PnbsfUkHAXQrRLSUlJbN++ncDAQAoLe/PWW30AWLjwIp6enbcdU0PCXQjR7pw7d47o\n6GiCg4NxcHBl/vz+lJVZc/fduYwfn6d1eRZBwl0I0a6cPn2aXbt2MXToUBwdHXnvvd6cOeNA795G\nFi26qHV5FkPCXQjRbhw7doz//ve/DB06FAcHB/bvd+Kjj9yxslJ47bVkHB2rtS7RYki4CyHahQMH\nDnD48GGGDx9Oly5dyM7W8/LL/VEUHX/4QwahoZ3vKtSmSLgLISyaoijs27ePEydOMHz4cOzs7Kio\n0PHiiwMoKLBh1KgCZs/O1rpMiyPhLoSwWIqisGfPHuLj4xkxYgR6vR6ANWs8OXXKEYOhnFdfTe70\nwx4bI+EuhLBI1dXVfPPNN6SkpDB8+HBsbNS4+vZbVzZsMGBtrbB8eRIuLlUaV2qZJNyFEBansrKS\n6OhosrKyuOWWW7C2tgYgPr4rS5d6A7BgQRrBwdJnvxEJdyGERamoqGDbtm3k5+czdOhQrH7uuVy5\nYsPzzw/EaLTivvuu8PDDlzWu1LJJuAshLIbRaOQ///kPZWVlhIaGotPpACgv1/GXvwwkJ8eWkJAi\nFi68yM9PiRto8jTE7NmzMRgMBNVbLjwvL4+IiAh8fX0ZO3Ys+fn5tc8tX76cQYMG4efnx7ffftt6\nVQshOpzS0lK2bNlCRUUFQUFBtcGuKPD66/04eVI9gfrmm4nY2so6zL+myXCfNWsWMTExDR6LjIwk\nIiKChIQExowZQ2RkJABnzpxh06ZNnDlzhpiYGObMmUN1tVxQIIT4dUVFRWzcuBErKysCAgIaPLdu\nnQfR0T2wt6/i7bcv0KNHpUZVti9Nhvsdd9yBq6trg8e2b9/OzJkzAZg5cyZbt24FYNu2bUybNg29\nXo+3tzc+Pj4cPny4lcoWQnQUBQUFbNy4ka5duzJ48OAGz33+uRvr1/fG2lohMjIJP79Sjapsf5o9\nOjQnJweDwQCAwWAgJycHgMzMTLy8vGpf5+XlRUZGhpnKFEJ0RHl5eWzYsAFXV1d8fHwaPPff/zqz\nYkVfAF56KZXf/OaaFiW2Wy06oarT6Wr7Yjd6vjFLliyp3Q4PDyc8PLwlZQgh2qHLly+zZcsWPD09\nGxwYAhw96siLLw6gulrHU09lct99uRpVqZ3Y2FhiY2NNfn+zw91gMJCdnY27uztZWVn06tULAE9P\nT9LS0mpfl56ejqenZ6OfUT/chRCdT1ZWFp999hn9+/end+/eDZ47dcqBBQt8MBqtuP/+yzz5ZJZG\nVWrrlwe+S5cubdb7m92WmTRpElFRUQBERUUxefLk2sc3btxIeXk5ycnJnD9/nhEjRjT344UQHVxa\nWhqbN2/Gx8fnumA/d86eZ5/1oaTEmnvuyWXRIhnyaKomj9ynTZvG3r17uXLlCn369OHVV19l0aJF\nTJ06lQ8++ABvb282b94MQEBAAFOnTiUgIAAbGxvWrl3bZMtGCNH5JCcns23bNgICAujRo0eD5y5c\n6MLcuYMoLLRh9OirLF6cws8XpgoT6BRFadMBozqdjjbepRDCApw7d46vv/6aoKAgXFxcGjx39qw9\nc+f61s7y+Pbbiej15smJyspK9u/fz4IFC8zyeVppbnbKFapCiFYXFxfHzp07CQ0NpVu3bg2eO3nS\ngWef9aGoyIbbb89nxYokswV7ZybhLoRoVcePH2fv3r21qyfVd+SIIwsW+FBaas1dd13ljTeSJdjN\nRMJdCNFqDh48yKFDhxg2bBj29vYNntu505VXXvGmosKKe+7JZfHiFGwkkcxGvkohhNkpisL333/P\n8ePHa1dPqnsOPv20F6tW9QHgoYcu8fzzaXLy1Mwk3IUQZnWj1ZMAKivhnXf6sGmTen3MvHnpPPpo\njgx3bAUS7kIIs6lZPSk5ObnB6kkA+fnWvPTSAA4fdkKvr2bJkhTGjbuqYbUdm4S7EMIsqqqqiI6O\nJjMzk2HDhtWungRw/rw9zz8/kIwMO7p3r2DFiiTCwoo0rLbjk3AXQrRYRUUF27dvJy8vr8HqSQDR\n0d1ZtqwvZWXWBAQU8+abibi7V2hYbecg4S6EaBGj0cjWrVspKSkhLCys9sr04mIr3nyzL199pV6J\neu+9ubz4YipdushQx7Yg4S6EMFlpaSlffPEFVVVVBAcH1z4eH9+Vv/61PxcvdsHOrpoXXrjIfffl\nyonTNiThLoQwSXFxMVu2bEGv1xMYGAioa52+/35v/v1vd6qqdAwaVMKyZcn071+mcbWdj4S7EKLZ\nrl27xubNm+nWrVvtIhunT3dl6VJvkpPt0ekUpk3L4ZlnMrCzkzaMFiTchRDNcvXqVTZv3oybmxve\n3t7k51uzdq0n//mPG4qio2/fMl55JYXQ0GKtS+3UJNyFWVRVQXo6ZGdDXh7k5qq3vDwoKYGKCvUC\nlpp7nQ7s7aFLl7qboyN07w49eqi37t3BYAAnJ63/60SNmtWTPDw88PDow5YtPVm3zoNr12ywtlb4\n3e+y+cMfMuWkqQWQcBfNkp8PJ06ot4QESEqCxERITVWDuzU4OSn06aOjTx/o0wf69lXvBwyAwYOh\nZ0/kRF0byM7O5rPPPqNfP2/i4/358589SE5W54sZMeIazz+fxoAB0lu3FBLu4oYUBeLjYd8+9Xbg\ngBrmN9K7N3h4NDzy7tFDPSK3sQEbm2pKSgopLMwlJ+cKV64UYWXVFSurruj1TpSX68nN1XH1qo6C\nAj1FRbYUFNhz7ZqeuDiIi2t8vy4u4Our3gYPVu+HDFHvZSIq80hPT+fzz7+goGAUa9f6Ex+vzu7o\n6Wlk3rx0Ro/Ol1+wFkb+1xcNXLsGu3bBV1/B119D1i+Wr+zSBYKCICwM/P1h4ED15u0NXbte/3l5\neXlkZGSQmJhIamoqjo42eHg48ZvfuNGjR4+fr2KsABpfAFlR1NZOUlIFqakKmZlWXLrUhcuX7bl0\nqRuXL/cgP9+Ww4fh8OGG77Wzg4AACA5Waw4OVm8Ggzm+qc7j/PkUliw5y4EDs0hOVudid3Mr54kn\nsrjvvlyZotdCyUpMguJi+PJL2LhRDfTy8rrneveG//kfuOMOuP12CAxs+mi4qKiIjIwMkpOTSUxM\npKqqCicnJ7p3707Pnj2xtbU1W92KolBUVEx6ejkXLliRmmpHeroDWVlOZGW5kZfn3Oj7evasC/rg\nYAgJUX8J1Ju4UACZmfDmm1eIirIlP1898eHmVs706Zd46KFL7aav3llXYpJw76QUBQ4ehH/+E7Zs\nUU96gtq7HjUK7r1XvQUFNd3PNhqNZGVlkZqaSmJiIgUFBTg5OeHi4kLPnj2vW5yhrRiNRrKzS4iP\ntyYhoQupqU6kpbmQmelGWdn1KW5jo+DnpyMkhAa3jnqUrygK5eXllJWVYTQaMRqNlJWVUVhYznff\n2bJ1qyv793enulqdRqB//1JmzMjh7rvzsLVtX/9+JdzbiIS7tkpKICoK1q2DU6fqHh81Ch55BKZM\nUY/Wb6SqqoqcnBzS0tJITEwkOzsbBwcHnJ2d6dmzJ87Ozha9MHplZRVJSZXExVmTkGBHcrIzFy+6\ncPmyC4pyfd29elUTGqojJKQu+AcPhnqz2GpKUZTacK4J6PrbJSUllJSUUFpaWrtd8xorKytsbGyo\nrrYjKakvR48O5Phxb0pL1b+urK2ruf32PB56KI8RIwqpN11MuyLh3kYk3LWRmwvvvQfvvgtXrqiP\n9eoFs2fDE0+offPGKIpCXl4eaWlpJCUlcfHiRWxtbXF2dsbNzQ1XV9cGs/+1V6WlVsTF6YiLs+Hc\nuS4kJ3cjLc2VsrLr20h6vYKfXxVhYVaEhVkREqL+hePmZvr+64f0LwO6JqTrB3RpaSmlpaWUl5dj\nbW2NjY0N1tbW6PV6rK2ta4Pb1tYWW1tb9Hp9vW1b0tO7ceyYEz/84MSPPzphNNYlt69vCWPH5jFh\nQi5ubpWm/0dZCAn3NiLh3rauXoWVK2H16rrWy/Dh8NxzcP/90FgLvLCwsEHfXFEUunXrRo8ePejZ\ns2eDxRc6supqyMy05exZO+LibEhIUEP/0qVujb7e2bkSH59qBg/WMXBgJf36GenTp4RevYqxti5r\ncCRdP6BLSkqoqKhoEM41gW1lZYVer68NZ71ej52dXYOff+0vpdxcG86d60pCgj2nTztw/Hg3Cgoa\nnjjx8yvmjjsKGDv2aoebKkDCvY1IuLeN4mI10N98EwoK1MfuvhsWLoQ772zYRzcajWRmZpKamsqF\nCxcoLCykW7duuLq60qtXL7o2NgymEysqsuLCBXsSEuyJj9eTkGBPaqojZWU3/qXn4FBGjx4luLmV\n0auXEYOhHHf3CtzclJ+HjSq4uFTh4FDdrCGF5eU6CgpsKCiwJj/fhuxsWzIy7EhPt6u9z8u7vq6e\nPcsJDS3i1luvMWrUNXr27LhT8HbWcDd5KOTy5cv55JNPsLKyIigoiA8//JDi4mIefvhhUlNT8fb2\nZvPmzbi4uJi6C2ECRYFNm+D55yEjQ31szBh44w0YOVL9uaqqiuzs7Nq+eU5ODo6Ojjg7OzNgwACc\nnRsfZSJUjo7VhIYWN7i8XlHUI+TU1C6kpnbh4kW72u3sbFuKi7tQXNyFixeb/mxra4UuXaqxs6vG\n1rYaW1sFW9tqFEVHZaWOqiqoqlK3CwutKSv79ZaYg0MVvr4l+PqW4u9fTGhoEZ6e5TIuvYMz6cg9\nJSWFu+66i/j4eOzs7Hj44YcZP348cXFxuLm58Ze//IUVK1Zw9epVIiMjG+5QjtxbzalT8Kc/wd69\n6s9Dh6pH7nfdpZCbm9ugb96lSxecnJxwc3Oje/fuDRZXEOalKHD1qg05ObZkZ9uSk6Ov3VaPutUj\n74ICG0pLm3f+wsamGmfnKpydK3FyqsRgqMDT04iXl7H2vmfPinZ7MtQc5Mi9GZycnNDr9ZSUlGBt\nbU1JSQkeHh4sX76cvT8ny8yZMwkPD78u3IX5GY3w2msQGanO8eLmBq+8Usro0SlcvJjMe+8l1fbN\n3dzcGDVqVIO1LUXr0umge/dKunevxN+/pMnXlpfrMBqtKC/XUV5uhdGo3ltZKVhbq0M2ra0VbGwU\nHByq6Nq1eW0c0XmY9C+8e/fuPPfcc/Tt2xd7e3vGjRtHREQEOTk5GH4eGGwwGMjJyTFrseJ6hw+r\nI17i4kCnU7jvvgzCw3dTWZnH4cPOuLq6Ehoair29vdalipugtmGqtC5DdAAmhXtiYiKrVq0iJSUF\nZ2dnHnroIT755JMGr9HpdDc8i79kyZLa7fDwcMLDw00po1OrrITXXlN4/XWortbh5naVWbP+yy23\nlNGrV1+cnIZoXaIQogViY2OJjY01+f0mhfuRI0cYNWoUPXqoayM+8MADHDhwAHd3d7Kzs3F3dycr\nK4tevXo1+v764S6aLy0Npk+Hfft06HQKDz6Ywrx5V+ja1Uvr0oQQZvLLA9+lS5c26/0mnWbx8/Pj\n4MGDlJaWoigKu3btIiAggIkTJxIVFQVAVFQUkydPNuXjRRO2b1evkty3D5ycivjb307x4ou5dO0q\njVchRB2TjtxDQkJ47LHHGDZsGFZWVgwdOpSnnnqKwsJCpk6dygcffFA7FFKYR3U1vPoq1PzyDgq6\nyMKFZ/Hz66FtYUIIiyQXMbUD167BjBnqUbuVFcyYEc/o0T8yZEig1qUJYfFkKKSwSMnJ6uyM8fHq\nohTLliVRWbmLgIBbtS5NCGHBOvGlDZbvyBG49VY12AMDISYmF6NxO6GhoXLRkRCiSZIQFio6Wp0D\n5tIldfqAPXuMnDr1H3x8fGSuFyHEr5Jwt0AffQSTJqmzOD72mBr0hw/vxNbWlt5NTbYuhBA/k3C3\nMGvXwqxZ6jQCf/2rGvRnz54kJSUFf39/rcsTQrQTEu4W5G9/g7lz1e233oLXX4crVy6zZ88egoOD\npc8uhLhpMlrGQixbph6pg7pi0pw56jzr27ZtY+DAgZqtRSqEaJ8k3C3A22+rwa7TwQcfqG0ZgJ07\nd2JjY4OHh4e2BQoh2h35O19j69apC2sArF9fF+ynTp0iOTmZgIAA7YoTQrRbEu4aiopS2y+gnkj9\n/e/V7cuX1T57SEhIh1h8WgjR9iTcNbJtmzoPO6gnT59+Wt0uLy9n+/bt9O/fX/rsQgiTSbhr4MAB\neOQRdTKwxYvhuefqntu1axfW1tZ4enpqV6AQot2TcG9jCQkwcSKUlcETT6jhXuP06dNcuHBB+uxC\niBaTcG9DOTlw992Qmwvjx6snU2sWq7py5Qq7d++WPrsQwiwk3NuI0QiTJ6uzPA4bBps2Qc0a1TV9\ndm9vbxwdHbUtVAjRIUi4twFFgT/+EQ4ehD59YMcOqJ/hu3fvxsrKCi8vWSZPCGEeEu5tYNUqdY6Y\nrl3VBTcMhrrn4uLipM8uhDA7CfdW9s03dRcpffQRhIbWPZebm8uuXbsIDg6WPrsQwqwk3FtRcnLd\nkMdXXoGHHqp7rqKiQvrsQohWI+HeSsrKYMoUyM9X52avP+QR1D67TqeTPrsQolVIuLeSBQvg2DHo\n31+dZqD+bL1nzpwhISFB+uxCiFYj4d4KPv0U/vEPsLODzz5TF7aukZuby86dO2U8uxCiVUm4m9nZ\ns/DUU+r26tUwdGjdc9JnF0K0FZPDPT8/nylTpuDv709AQACHDh0iLy+PiIgIfH19GTt2LPn5+eas\n1eIZjTBtmrr26fTpdSFfY8+ePQDSZxdCtDqTw33evHmMHz+e+Ph4Tp48iZ+fH5GRkURERJCQkMCY\nMWOIjIw0Z60W7//9PzhxAgYMUKfwrZlaACA+Pp6zZ88SGBioXYFCiE5DpyiK0tw3FRQUEBYWRlJS\nUoPH/fz82Lt3LwaDgezsbMLDwzl79mzDHep0mLBLi7drF0REgLU1fP893Hpr3XN5eXl8/PHHhIaG\nSjtGiDZWWVnJ/v37WbBggdaltEhzs9OkI/fk5GR69uzJrFmzGDp0KE8++STFxcXk5ORg+PnyS4PB\nQE5Ojikf3+7k5sLMmer24sUNg72mz963b18JdiFEmzFpDdXKykqOHTvG3//+d4YPH878+fOva8Ho\ndDp09fsS9SxZsqR2Ozw8nPDwcFPKsBhPPw2ZmfCb38CLLzZ8LjY2lurqavr27atNcUKIdik2NpbY\n2FiT329SWyY7O5vbbruN5ORkAL7//nuWL19OUlIS3333He7u7mRlZTF69OgO35bZsgWmTlUnAjt5\nUh3XXuPs2bN8++23jBw5EhsbWYtcCC1IW6YZ3N3d6dOnDwkJCYC6elBgYCATJ04kKioKgKioKCZP\nnmzKx7cbly7VrYH61lsNgz0vL49vvvmG4OBgCXYhRJszOXXeffddpk+fTnl5OQMHDuTDDz+kqqqK\nqVOn8sEHH+Dt7c3mzZvNWavFeeYZuHIFxoxpOOyxfp+9W7du2hUohOi0TA73kJAQfvzxx+se37Vr\nV4sKai+2bFFvjo7wr381HPa4d+9eqqqqpM8uhNCMXKFqgtxcmDtX3V65Ery96547e/YscXFxDBky\nRJPahBACJNxN8pe/wOXLMHo0/OEPdY9fvXpV+uxCCIsg4d5Me/fC+vVga6tODlbTjqmsrOTLL7+k\nT58+ODk5aVukEKLTk3BvBqNRXQsV4KWXwNe37rm9e/dSUVFBv379tClOCCHqkXBvhpUr1VkffX1h\n0aK6x8+dO8fp06elzy6EsBgS7jfp/Hl4/XV1u2audpA+uxDCMkm436T589W2zGOPqSdSQe2z79ix\nA09PT+mzCyEsioT7TfjqK4iOBicntTVT47///S9GoxHv+mMhhRDCAki4/wqjUT1qB1i6FHr1UrcT\nEhI4deoUQUFB2hUnhBA3IOH+K1atggsXwN+/7sKl/Px8YmJiCAoKkj67EMIiSbg3ITMTXntN3V69\nGvR6qKqq4ssvv8TT0xNnZ2dtC+xkNm/ezJ133snp06e1LkUIiyfh3oSFC6G4GO6/X11lCWDfvn2U\nlZVJn10D9957L3Z2drJUoRA3QcL9Bo4cgU8+Ua9Efftt9bELFy7w008/ERwcrG1xndSRI0cICwu7\n4SIwQog6Eu6NUBR44QV1e948dZ72goICoqOjGTJkiPTZNXLw4EF0Oh0xMTEsW7aMCxcuaF2SEBZL\nwr0RO3ZAbCx0765OM1BVVVU7nt3FxUXr8jqFjRs3MmbMGB599FFSU1MBOHz4MNOnT+fuu+/mf/7n\nf1i7dq3GVQphuSTcf6GyUp31EeDll8HFRe2zl5SUSJ+9jRw5coR33nmHVatWUVxczGuvvUZ2djaK\notQOPb0OCbJuAAAVJklEQVRy5Qr5+fkaVyqE5ZJw/4UPPlDnjxkwQF1CLzExUfrsbWzNmjXcdttt\n+Pr6oigKBoOB+Ph4QkNDa19z6NAhRo0apWGVQlg2Cfd6iorglVfU7chIKC0t4KuvvpI+exs6ffo0\n8fHxREREYGdnx9atW3njjTdwcHCoXbLw4sWLXLhwgUcffVTjaoWwXBLu9axapS56PXIk3H+/2mf3\n8PCQPnsb+uqrrwCuOyofPnw4VlZW7Nixgw0bNrBu3Tq6dOmiRYlCtAs6RVGUNt2hTkcb7/KmXL1a\nMyoGdu8GK6tYzp8/z9ChQ7UurVO59957cXR0ZNOmTVqXIjqIyspK9u/fz4IFC7QupUWam51y5P6z\nlSvVYB8zBry9kzhx4oTMG9PGLl68yKVLlxr01oUQppFwB3Jy1OkFAF56qZjo6GiCgoLQ6/XaFtbJ\n/PjjjwCy6IkQZiDhDixbBiUlMGGCQlbWVgwGg/TZNXDkyBEA/P39Na5EiPavReFeVVVFWFgYEydO\nBCAvL4+IiAh8fX0ZO3ZsuxiHfPGiurISwAMPHKOkpIQBAwZoW1QndezYMWxtbenfv7/WpQjR7rUo\n3FevXk1AQEDtXB+RkZFERESQkJDAmDFjiIyMNEuRrWnZMigvhwkTiigo2Cd9do2kpqaSl5fHwIED\nsba21rocIdo9k8M9PT2d6OhonnjiidozuNu3b2fmzJkAzJw5k61bt5qnylaSlgbr14NOpxASspUh\nQ4ZIn10jx48fB2Dw4MEaVyJEx2ByuC9YsICVK1diZVX3ETk5ORgMBgAMBgM5OTktr7AVrVgBFRVw\n222pBAfb4urqqnVJndaxY8cA8PHx0bgSIToGky673LFjB7169SIsLIzY2NhGX6PT6W44NeuSJUtq\nt8PDwwkPDzeljBbJyID331e3IyJ+ZODAgW1eg6hz6tQpwDLCvaqqyuTWUGVlpVzNLMwiNjb2hvl6\nM0z6v3D//v1s376d6OhoysrKuHbtGjNmzMBgMJCdnY27uztZWVn0qllw9Bfqh7tWVq5Ue+0hIee5\n++4+WpfTqV29epX09HR0Op3mv2T37NlDcXFx7SCB5vrwww8ZMWIEISEhZq5MdDa/PPBdunRps95v\nUltm2bJlpKWlkZyczMaNG7nrrrv4+OOPmTRpElFRUQBERUUxefJkUz6+1WVnwz//qZ4nePrpXOmz\na+zkyZMAuLq6tskQ1LS0NObPn8+aNWtYvnx57Tmjo0ePcvz4cZODHWDWrFmsX7+e5OTkm37PO++8\nw7333svw4cM5evSoyfsWoj6zjHOvab8sWrSInTt34uvry549e1i0aJE5Pt7sVq5UKCvTMXx4OsOG\nSbBrrSbc26IlU1FRwTPPPMOYMWPIzc1l27ZtFBcXU1RUxJo1a3jmmWda9Pk2Nja8+OKLLF68mMrK\nypt6z4IFC5g5cya2trYyWkuYTYubg3feeSd33nknAN27d2fXrl0tLqo15ebCunUKoGPevGtalyOo\n67cPGjSo1fd14MABMjMzGTp0KAMGDKidy+bdd9/lnnvuwc7OrsX7cHd3Z+DAgezYseOm/3o9fvw4\nAQEB2Nratnj/QkAnvEL1vfegtNSKoKB0/PxKtS6n06uqqiI+Ph5omyP3o0eP4urqiqenJ4GBgYwY\nMYLS0lK2bt3K+PHjzbafhx9+uLZFeTNOnDghk9QJs+pU4V5cDGvWqNvjx5/WthgBQEpKCmVlZeh0\nOnx9fVt9f3FxcQQEBDR47Pvvv8fDwwMnJyez7cfX15eCggLOnj37q69NT0/nypUrEu7CrDrVmK31\n69W2TFBQKYMHZwM9tS6p0ztz5gwA1tbWrTrtw5IlS8jLy+Onn37C29ubZ599Fk9PTxYuXMihQ4ea\nXGkrPj6e6OhorKysyMzM5OWXX+aLL76gsLCQS5cu8Yc//AEvL68G77GysiIkJISDBw/i5+fX4Lkf\nf/yRL774Ag8PDwoLC2uvyv3lCBtT9itEjU4T7hUV8Pbb6vYTT+RygyH4oo3VhHv//v1bdXz4kiVL\nyMjIYPLkycydO7fBELOEhATuv//+Rt+Xnp7O9u3bWbhwYe3n/P73v2fp0qVUV1fz5JNP4ufnx/Tp\n0697b79+/UhISGjw2NatW1m3bh2ffPIJPXv2JDs7mwcffJCAgIAGi4+0ZL9CQCdqy2zeDKmp4OsL\nd91VqHU54mc14f7Lo9vWcO7cOYDr2j+ZmZm1S/j90qeffsqzzz5b+3NpaSnOzs4EBQXh7u7O9OnT\nbzh0slu3bmRmZtb+nJCQQGRkJM899xw9e6p/Nbq7u2Nvb88tt9xitv0KAZ0k3BUF3nxT3f7LX0Dm\npbIMVVVVXLhwAWibaX4TEhJwdHTEw8OjweNFRUU4Ojo2+p4ZM2Zgb29f+/OpU6cYMWIEoE6xMW/e\nvBv26p2dnSkqKqr9ee3atTg4ODBmzJjax5KSkigoKLiu396S/QoBnSTcv/0WTp4EDw+QNZUtR0pK\nCuXl5eh0ujYJ93PnzjV60rap5cvq/yJISUnh8uXLDBs27Kb2pyhK7ecWFhZy4MABRo4c2WBqg6NH\nj2JlZXXd6lMt2a8Q0EnC/W9/U+//9CcwwzBmYSY1/WgbG5s2mQ3y/PnzjYZ7t27dKCgo+NX3Hzly\nBL1e3+Dka3p6+g1fX1BQUNvuSUtLo7q6+roTt0eOHMHf3x97e3syMjLMsl8hoBOEe1yceuTetSs8\n9ZTW1Yj6asLdx8en1aeAKCgoICcnp9Fw9/DwaDTcy8rKWL16dW3r6NChQwwaNKj2Qqfq6mo+/vjj\nG+7z2rVreHp6AuDg4ACoPfb6n3/s2LHalszGjRvNsl8hoBOMllm1Sr2fORO6d9e2FtFQYmIiAIGB\nga2+r5qTqY1dBRsaGtroXDA//PADn3zyCf7+/tjY2JCWltbgxOv69eubPKmZnJzMyJEjAejbty+D\nBg2qPTqvrKxkxYoVVFRU4OXlRV5eHt1//h+0pfsVAjp4uF++DDUHOPPmaVuLuF7NkWlbhPvZs2fp\n1q1bo0fut912G2/XjJOt55ZbbmHChAnEx8dz7tw5PvroIyIjI1m2bBl6vZ4777zzhot5V1ZWcvLk\nydoRLzqdjsjISP72t7+Rk5NDdXU1jz/+OLfccgs7duwgPj6+9rUt2a8QNTp0uP/jH2A0wr33gizw\nY1kKCwu5fPkyOp2uTYLq7NmzDB8+vMHiMjXCwsLIzc3l8uXLtUMUAVxcXFi8eHGD197sdNVxcXEY\nDIYGfyn07duXVTV/Sv7My8uLCRMmNHisJfsVokaH7bkbjeo8MgALFmhbi7heTUvGyckJb2/vVtnH\nRx99xNy5cwF1PH39IYj12draMnXqVDZs2GC2ff/f//0fj8rQLKGhDhvuGzZATg4EB8Ndd2ldjfil\npKQkgOuGAJrT119/ja2tLefPn0ev198w3EFd83f//v1cu9bymUJTUlLIzs6WvrjQVIcMd0WBd99V\nt+fNQ6YasEA14R4WFtZq+5gxYwZubm6sX7+elStXNrl0XpcuXXj55Zd5/fXXbzjm/WYYjUZWrlzJ\nG2+8ccNlJoVoCx2y537oEBw7Bj16wLRpWlcjGlMzDLI1j9wnTJhwXT+7KYGBgTzwwANs2rSJRx55\nxKR9fvjhh8ydO1cm9BKa65Dh/ve/q/ePPw71ruAWFuTChQvY29u3yZwyzXHrrbdy6623mvz+P/7x\nj2asRgjTdbi2zKVLsGWL2op5+mmtqxGNycrKorCwkCFDhjTZKhFCmK7Dhfu//gXl5TBhArTSIAzR\nQjUrL/1yJkQhhPl0qHCvrIR169Ttn0fACQsUFxcHUDvLoRDC/DpUuH/5JaSnw6BBEBGhdTXiRk6d\nOoWDg0ObXJkqRGfVocK95qKlOXOgkQsRhQUoKysjLi6OkSNHNnq1qBDCPDrMv66EBNi9W5398fe/\n17oacSNHjhyhvLycO++8U+tShOjQTAr3tLQ0Ro8eTWBgIEOGDGHNmjUA5OXlERERga+vL2PHjiU/\nP9+sxTblX/9S7x95BFxc2my34le89dZbTJs2jcrKSgBiYmJwcnJq8mpRIUTLmRTuer2ed955h7i4\nOA4ePMh7771HfHw8kZGRREREkJCQwJgxY4iMjDR3vY0yGuHDD9VtmbPdshw+fJiysjKqq6vJzs5m\nz549TJs2rXZuciFE6zDpIiZ3d/faRQccHR3x9/cnIyOD7du3s3fvXkCdqyM8PLxNAn7bNrhyRZ1H\nRgZgWJaQkBC6d+/OtWvXePXVV+nbty+/l76ZEK2uxT33lJQUjh8/zsiRI8nJycFgMADqIr45OTkt\nLvBm/POf6v1TT8k8MpZm7ty5xMXFMXnyZGxtbXn33XexsWn8mKKyspJ169bx2WefsXHjRhYsWCDL\nyQlhohZNP1BUVMSDDz7I6tWrG6wUA+riBDeaOKn+3NTh4eGEh4ebXMP587BnjzrNwPTpJn+MaCUu\nLi78vWY+iF+xfPlyBg0axJQpU8jPz+d///d/ZY4W0WnFxsYSGxtr8vtNDveKigoefPBBZsyYweTJ\nkwH1aD07Oxt3d3eysrLo1atXo+8158IDNSdSH35YTqS2Z+fPn2fnzp0sXLgQUOeeqVlbVIjO6JcH\nvkuXLm3W+01qyyiKwuOPP05AQADz58+vfXzSpElERUUBEBUVVRv6raW8XE6kdhSHDh0iNDQUW1tb\nQD0RO2zYMAoLCzWuTIj2yaRwr1nA97vvviMsLIywsDBiYmJYtGgRO3fuxNfXlz179rBo0SJz19vA\n9u3qOqlDhkALJvITFsDZ2ZkePXoAUFJSwnfffUdISAhff/21xpUJ0T6Z1Ja5/fbbqa6ubvS5Xbt2\ntaig5li/Xr1/4gk5kdrejRs3jhMnTvDNN99QXl7OuHHj+OGHHyxuSmAh2ot2O597RgZ88w3o9XIi\ntSOwtbXl5Zdf1roMITqMdjv9wMcfQ3U1TJoEbm5aVyOEEJalXYa7otSdSJ01S9tahBDCErXLcD9w\nQJ0orHdvGDdO62qEEMLytMtwrzmROmMG3OBiRyGE6NTaXbgXF8OmTeq2tGSEEKJx7S7cP/8ciorU\nce0ySk4IIRrX7sL95wtg5ahdCCGa0K7CPT0dvvsO7OzUuWSEEEI0rl2F+4YN6jDIiRPB2VnraoQQ\nwnK1q3D/9FP1Xq5IFUKIprWbcD99Gn76CVxd4Z57tK5GCCEsW7sJ95qj9oceUnvuQgghbqxdhHt1\ndV24P/qotrUIIUR70C7C/fvvIS0N+vaF3/xG62qEEMLytYtw/+QT9X76dLBqFxULIYS2LD4qjUbY\nskXdllEyQghxcyw+3KOjIT8fQkMhMFDraoQQon2w+HCXse1CCNF8Fh3u+fnw5Zfq+qjTpmldjRBC\ntB8WHe6ffw7l5TB6NHh6al2NEEK0HxYd7jK2XQghTGP2cI+JicHPz49BgwaxYsUKkz8nJwf27gVb\nW7j/fjMWKIQQnYBZw72qqopnnnmGmJgYzpw5w4YNG4iPjzfps774Qr0ydexYcHExZ5WW48iRI1qX\nYDHku6gj30Ud+S5MZ9ZwP3z4MD4+Pnh7e6PX63nkkUfYtm2bSZ9VM7Z96lQzFmhhjh49qnUJFkO+\nizryXdSR78J0Zg33jIwM+vTpU/uzl5cXGRkZzf6c+i2ZSZPMWaEQQnQONub8MJ1OZ5bPqd+SaY1F\nOXQ6HXl5eRw7dsz8H94MWVlZmtdgKeS7qCPfRR1zfBfV1dXY2Jg16toHxYwOHDigjBs3rvbnZcuW\nKZGRkQ1eM3DgQAWQm9zkJje5NeM2cODAZuWxTlEUBTOprKxk8ODB7N69Gw8PD0aMGMGGDRvw9/c3\n1y6EEELcBLP+rWJjY8Pf//53xo0bR1VVFY8//rgEuxBCaMCsR+5CCCEsQ5teoWquC5zau7S0NEaP\nHk1gYCBDhgxhzZo1WpekuaqqKsLCwpg4caLWpWgqPz+fKVOm4O/vT0BAAAcPHtS6JM0sX76cwMBA\ngoKC+N3vfofRaNS6pDYze/ZsDAYDQUFBtY/l5eURERGBr68vY8eOJT8/v8nPaLNwN+cFTu2dXq/n\nnXfeIS4ujoMHD/Lee+912u+ixurVqwkICDDbiKv2at68eYwfP574+HhOnjzZaduaKSkpvP/++xw7\ndoxTp05RVVXFxo0btS6rzcyaNYuYmJgGj0VGRhIREUFCQgJjxowhMjKyyc9os3A35wVO7Z27uzuh\noaEAODo64u/vT2ZmpsZVaSc9PZ3o6GieeOIJOnOXsKCggH379jF79mxAPYfl3BpjgdsBJycn9Ho9\nJSUlVFZWUlJSgmcnmj3wjjvuwNXVtcFj27dvZ+bMmQDMnDmTrVu3NvkZbRbu5rrAqaNJSUnh+PHj\njBw5UutSNLNgwQJWrlyJVSdfQzE5OZmePXsya9Yshg4dypNPPklJSYnWZWmie/fuPPfcc/Tt2xcP\nDw9cXFz47W9/q3VZmsrJycFgMABgMBjIyclp8vVt9q+ps/+53ZiioiKmTJnC6tWrcXR01LocTezY\nsYNevXoRFhbWqY/aQR1KfOzYMebMmcOxY8dwcHD41T+9O6rExERWrVpFSkoKmZmZFBUV8WnNNLEC\nnU73q5naZuHu6elJWlpa7c9paWl4eXm11e4tTkVFBQ8++CCPPvookydP1roczezfv5/t27fTv39/\npk2bxp49e3jssce0LksTXl5eeHl5MXz4cACmTJnSaa9UPXLkCKNGjaJHjx7Y2NjwwAMPsH//fq3L\n0pTBYCA7OxtQr9zt1atXk69vs3AfNmwY58+fJyUlhfLycjZt2sSkTjpxjKIoPP744wQEBDB//nyt\ny9HUsmXLSEtLIzk5mY0bN3LXXXfx73//W+uyNOHu7k6fPn1ISEgAYNeuXQR20oWD/fz8OHjwIKWl\npSiKwq5duwgICNC6LE1NmjSJqKgoAKKion79oLBF8w00U3R0tOLr66sMHDhQWbZsWVvu2qLs27dP\n0el0SkhIiBIaGqqEhoYqX3/9tdZlaS42NlaZOHGi1mVo6sSJE8qwYcOU4OBg5f7771fy8/O1Lkkz\nK1asUAICApQhQ4Yojz32mFJeXq51SW3mkUceUXr37q3o9XrFy8tLWb9+vZKbm6uMGTNGGTRokBIR\nEaFcvXq1yc+Qi5iEEKID6tzDE4QQooOScBdCiA5Iwl0IITogCXchhOiAJNyFEKIDknAXQogOSMJd\nCCE6IAl3IYTogP4/oENi/a9Z/A8AAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute the integral both at high accuracy and with the trapezoid approximation" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from __future__ import print_function\n", "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\n" ] } ], "prompt_number": 5 } ], "metadata": {} } ] }