{
"metadata": {
"name": ""
},
"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": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAF9CAYAAADY7xbEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8XXWd//HXufdmz02atEma7mlLS6F0BdlK2RcR0NHK\npugwKi7g/H6MMw/8OYAOHWVwm2GUEUdGHcsiAq6IVNoqFESWphvd9zT3JrnZ776dc35/3CY0NE3b\n5CY3N/f9fDzySHLu9977KYnt2+/5nM8xbNu2EREREZG0cWS6ABEREZGxRgFLREREJM0UsERERETS\nTAFLREREJM0UsERERETSTAFLREREJM0UsERERETSzDXQg8lkkieffJIdO3bw4IMP9h63LIuHH34Y\n27b5h3/4BwC2bNnCM888g9PppKKigjvvvBOXa8CXFxERERmTBtzBeuqpp5g/f/4xx5977jnOO++8\n3u9jsRhPPPEE9957L1/72tdYvHgxzz33XPqrFREREckCAwas2267jSVLlvQ59uabb1JeXs7s2bN7\nj+3YsYMlS5ZQUFAAwLJly3jnnXeGoVwRERGR0e+UerC8Xi+bN2/mqquu4ug77ASDQdxu97sv6nBg\nWVb6qhQRERHJIqfUJLV161ZaW1v51re+RTwe5/Dhw6xatYqzzjoLj8fTu86yLByOgbPb2rVrB1ex\niIiISAZcfvnlJ732lALW1VdfzdVXXw1Aa2srq1at4rbbbiMajfLzn/+cD37wgxQWFvLqq6/227v1\nXu89/SjZ43e/+x3XX399psuQQdDPLrvp55e9cvVnd7Azwh3P7STPabDqpjOpLM7LdEmDUl9ff0rr\nB32Zn23bGIYBQGFhIbfeeisrV67E5XJRUVHBXXfdNdiXFhERkTHi6c0tAFwzZ3zWhqvBOKmAdfSI\nhh7V1dXcfffdvd8vWLCABQsWpK8yERERyWpN/hh/2teJw4CPLqjOdDkjSoNGRUREZFg8s9WHZcPl\nsyuZ6C7IdDkjSgFLRERE0q49nGD1rnYM4KYFNZkuZ8QpYMmgzJkzJ9MlyCDpZ5fd9PPLXrn2s3tm\nSwsJy+bCGeVMqyjMdDkjTgFLBmXu3LmZLkEGST+77KafX/bKpZ9dRzjB8zvaALh10cQMV5MZClgi\nIiKSVr/Y0kLctLlgejmzJxRnupyMUMASERGRtGk/avfqtiW5uXsFClgiIiKSRr/YnNq9unBGObPG\n5+buFShgiYiISJq0hxP8fmdq9+rji2szXE1mKWCJiIhIWjx9ZPdq2YxyZo0vynQ5GaWAJSIiIkPW\nFoq/u3u1JLd3r0ABS0RERNLg6c0+EqbNRXXjmFmZ27tXoIAlIiIiQ9QWivPCrp7eq9y9cvBoClgi\nIiIyJE9vbiFh2iyvG0eddq8ABSwREREZAl8wzgs7U/cc1O7VuxSwREREZNBW1TeRsGwunlnBDO1e\n9VLAEhERkUFp6Izy0p4OnAZ8cqmuHDyaApaIiIgMyk83eLFseP/pE5hcXpDpckYVBSwRERE5ZTt9\nIV492E2B0+Bji9R79V4KWCIiInJKbNvmf97yAvCh+dWML8nLcEWjjwKWiIiInJJ6b4DNTUFK853c\nuKA60+WMSgpYIiIictJs2+bHR3avblxYg7vAleGKRicFLBERETlp6w92sactQmWxiw+dWZXpckYt\nBSwRERE5KaZl89O3mwD4+OJaCl2KEcej/zIiIiJyUv64u53G7hiTyvK5Zu74TJczqilgiYiIyAlF\nEib/W5/avfrk0lpcDiPDFY1uClgiIiJyQs9s8dERTjK3qpiLZ1ZkupxRTwFLREREBtQaivPMlhYA\nPnvuZByGdq9ORAFLREREBvTTt5uImTYX1Y1j/sTSTJeTFRSwRERE5Lh2t4V5aU8HeQ6DT50zKdPl\nZA0FLBEREemXbdv89xseAD50ZhWTynRD55OlgCUiIiL9+suhbrY0BSkrcHLLoppMl5NVFLBERETk\nGAnT4rE3U7fEuW1JLaW6Jc4pUcASERGRYzy/ow2PP8bU8gI+MG9CpsvJOgpYIiIi0oc/muTxjc0A\nfObcyRoqOggKWCIiItLH/25oIhAzWTSplHOnlmW6nKykgCUiIiK99raF+f3ONhwGfOG8KRgaKjoo\nClgiIiICgGXbfP/1Riw7NZZhRmVRpkvKWgpYIiIiAsDavZ1sbwlRUeTitiW1mS4nqylgiYiICKG4\nyWNvpoaKfuZ9kynJd2a4ouymgCUiIiKsqm+iM5LkzJoSLp9dkelysp4CloiISI470BHh19tacRhw\n1wVqbE+HAceyJpNJnnzySXbs2MGDDz4IwGuvvcYLL7xAXl4eJSUlfOELX6CkpIQtW7bwzDPP4HQ6\nqaio4M4778Tl0tRXERGR0cy2bR450th+wxkTmDW+ONMljQkDJqCnnnqK+fPns2PHDgBM06S+vp6v\nfvWr5Ofns3r1atasWcM111zDE088wQMPPEBBQQGvvPIKzz33HDfddNOI/CFERERkcP68v5MtTUHK\nC118cqka29NlwFOEt912G0uWLOn93ul08sUvfpH8/HwAotEoVVVV7NixgyVLllBQkLrL9rJly3jn\nnXeGsWwREREZqlDc5IdvpBrbP3XOJNy632DaDLoH65VXXqGxsZHzzz+fYDCI2+1+90UdDizLSkuB\nIiIiMjx+/JaXjnCSedXFXDWnMtPljCmDiqrPPfcc4XCYO++8E4DS0lI8Hk/v45Zl4XCof15ERGS0\n2tYS5Hc72nAacPeyaTjU2J5Wp5SCTNPkhz/8ISUlJdx22229x+fNm8fGjRuJRqMAvPrqq8yfPz+9\nlYqIiEhaJEyL/3j1MAA3LazRxPZhcEo7WBs3buTNN9+kubmZN954A4AlS5Zw/fXXc+utt7Jy5Upc\nLhcVFRXcddddw1KwiIiIDM0zW3wc6owyuayAWxdNzHQ5Y5Jh27adiTdeu3ZtnwZ6ERERGX6N3VE+\n+8udJEybb147m0WT3Cd+klBfX8/ll19+0uvVKCUiIpIjbNvm4VcPkzBtrjqtUuFqGClgiYiI5Ig/\n7ulg85GZV3ecOznT5YxpClgiIiI5oCuS4L+PzLz63HmTKSvUzKvhpIAlIiKSAx75SyOBmMnSyW4u\nm6WbOQ83BSwREZEx7uX9nbx8oIuiPAf/Z9lU3cx5BChgiYiIjGGdkQTfey018+qO901morsgwxXl\nBgUsERGRMarnqkF/zGTJZDfXnj4+0yXlDAUsERGRMWrdvk7+cqib4jwH/3DRNJ0aHEEKWCIiImNQ\neyjBI39JnRr83HlTqC7Nz3BFuUUBS0REZIyxbZtvv3yAYNxiQVU+V8+pzHRJOUcBS0REZIz57RYP\nG7whipzwdwvKdWowAxSwRERExpBdjT7+p74VgL+Z4aSi0JnhinKTxriKiIiMAbZt09zi47uvNRE1\nYX6FwdIJ2rnKFAUsERGRLGeaJs3NzTy9tY0DAZuyPLhxplOnBjNIAUtERCSLJRIJvF4vW7zdrG60\nMIBbZzspzVO4yiQFLBERkSwVi8VobGykIxDmib0mNnDZJAdzytVinWkKWCIiIlkoHA7T2NhILBbj\n2QMmHTGYUmJwzRSFq9FAAUtERCTL+P1+PB4Ppmmyoc1mY7tNvgM+PtuJy6FTg6OBApaIiEgW6ezs\nxOv1Yts2bVGb5w6YAHx4hpPqIoWr0UIBS0REJAvYtk1bWxstLS0AJC2bVXtMYhYsrDQ4p0rhajRR\nwBIRERnlTNOkpaWFjo6O3mO/PmRxOGRTWQAf1UiGUUcBS0REZBRLJBI0NTXh9/t7j21os/hLi4XT\ngE+e5qLYpXA12ihgiYiIjFKxWAyPx0M4HO491hy2eWb/u31XU0sVrkYjBSwREZFRKBKJ9I5h6BE1\nbX66O0ncgrMnGJxXrXA1WilgiYiIjDKBQACPx0Mymew9Zts2v9hv4ovCxCJYUae+q9FMAUtERGQU\nOXoMw9FebbHY1G5T4IS/neMi36lwNZopYImIiIwCPWMYfD7fMeHqQMDit4csAG6eqXlX2UABS0RE\nJMMsy6KlpYX29vZjHuuM2fxkt4lpw8UTHSwcr1vhZAMFLBERkQxKJpM0NTXR3d19zGMx0+bHu5IE\nEzCn3OC66QpX2UIBS0REJEPi8Tgej4dQKHTMY7Zt8/N9Jp4wTCiET5zmxKmm9qyhgCUiIpIB/Y1h\nONpLHovNHTaFTvjUXA0TzTYKWCIiIiMsFApx+PDhPmMYjra53eLFRgsDuO00JzVqas86ClgiIiIj\nqLu7G4/Hg2VZ/T7uCdk8tS81qf36aQ7mjVPfVTZSwBIRERkhbW1ttLS0HDOGoUd33OZ/dr07qf3i\nWoWrbKWAJSIiMswsy8Ln89HW1nbcNdGkzY92JumKw4xSg4/O1KT2bKaAJSIiMoySySTNzc10dXUd\nf42VmnXlDUNVIXxqrpM8h8JVNlPAEhERGSbxeByv10swGDzuGtu2eXq/yR6/jTsP7jjdRUmewlW2\nU8ASEREZBpFIBI/HQzQaHXDdC4ctNrTZ5Dvg03NdjC9UuBoLFLBERETSLBQK0djYSCKRGHDday0m\na70WDuCTc5xMLVW4GisUsERERNLI7/fj8XgwTXPAdVs7LH55IDWq4aMznRrHMMYoYImIiKRJR0cH\nTU1Nxx3D0GNnl8XP9pjYwNVTHJxbrXA11ihgiYiIDJFt27S2tuLz+U64dp/f4se7TEwblk90cNVk\nhauxaMCAlUwmefLJJ9mxYwcPPvggAOvXr2f16tU4HA7q6uq4/fbbAdiyZQvPPPMMTqeTiooK7rzz\nTlwu5TcRERnbTNOkubmZzs7OE649FLB4bKdJ0obzqg0+ON2hWVdj1ICx+amnnmL+/Pm93/t8Ptau\nXcvKlSt54IEHcLvdrFu3jlgsxhNPPMG9997L1772NRYvXsxzzz037MWLiIhkUiKRoLGx8aTClSdk\n8987TWIWLBlvsKJOg0THsgED1m233caSJUt6v9+0aRPLly/v/YW44oor2LBhAzt27GDJkiUUFBQA\nsGzZMt55551hLFtERCSzotEoDQ0NBAKBE65tidj8cEeSiAnzKwxumeXEoXA1pp3Sid9gMEhZWVnv\n92VlZfj9foLBIG63+90XdTiOexNLERGRbBcKhTh06BCRSOSEa30Rmx9sTxJMwtxyg0+c5sSpKe1j\n3ikFLLfbjd/v7/3e7/dTVlZGaWlpnwRvWRYOh5r2RERk7PH7/TQ0NJxwxhVAc9jmke1J/AmY6Ta4\nfY4Tl8JVTjilFLRo0SLWr1/fuzu1Zs0ali5dyrx589i4cWPvtNpXX321T++WiIjIWNDZ2cnhw4dP\nOOMKUj1Xj2xPEkjAaWUGnzndSb5T4SpXnNJlflVVVVx66aXcd999OJ1Opk+fzooVKwC49dZbWbly\nJS6Xi4qKCu66665hKVhERGSk2bZNW1sbLS0tJ7W+IWjxwx0mERPmjTP42zm6eXOuMewTTUMbJmvX\nru3TQC8iIjIamaZJS0sLHR0dJ7X+QMBKXS14pKH9E6dl9rTg+PHjqa2tzdj7jxX19fVcfvnlJ71e\ng6pERESOI5FI0NTU1Kf/eCB7ui3+Z5dJ3IKFlQYfn62G9lylgCUiItKPWCyGx+MhHA6f1PpN7RZP\n7E1NaD97gsFNs5w4NYohZylgiYiIvEckEqGxsZFYLHZS619pMvnNIQsbuLDGwd/McGjOVY5TwBIR\nETlKIBCgsbHxpK4UtGyb5xss/tyUurr+A1MdXDZJt78RBSwREZFenZ2deL1eTub6r6Rl8/N9JvXt\nNg4Dbp7p5OwqzYCUFAUsERHJeac6hiGStPnpbpM9fpsCB/ztHCdzxylcybsUsEREJKdZlkVLSwvt\n7e0ntd4XsfnxriS+KLjz4DOnu5hSolOC0pcCloiI5KxkMklTUxPd3d0ntX5nl8XP9phETagthk/N\ndVFZoHAlx1LAEhGRnBSPx/F4PIRCoROutW2bPzdZPN+QulJwQaXBLbOcFOjWN3IcClgiIpJzTmUM\nQ8Ky+cV+kw1tqcb3q6c4uHKyxjDIwBSwREQkpwSDQRobG0kmkydc2x61+dkek8Mhm3wH3DrbyYJK\nNbPLiSlgiYhIzujq6sLr9WJZ1gnXbm63eHp/qt+qsgD+bo6LSWpml5OkgCUiImOebdu0t7fT0tJy\nwhlXCcvmt4csXmtJhbD5FQY3z3JS7FK4kpOngCUiImOaZVn4fD7a2tpOuLY1YvOzPUk8YXAacP00\nBxdN1GR2OXUKWCIiMmYlk0mam5vp6uoacJ1t27zdZvPLAyYxC8YXwCdOczK1VP1WMjgKWCIiMibF\n43G8Xi/BYHDAdf64zTMHTLZ1pk4dLqw0uHGmkyKdEpQhUMASEZExJxKJ4PF4iEajA67b2Gbxy4Mm\noSQUOuFD052cU2XolKAMmQKWiIiMKSczhiGYsPnlQZNN7aldqznlBjfNdFKhqeySJgpYIiIyZnR3\nd+P1ejFNs9/HbdtmQ5vNbxtMggnId8AN0x2cX61GdkkvBSwRERkT2tvbaW5uPu4YhuawzXMHTfb5\nU4/PKjO4eaaT8YUKVpJ+ClgiIpLVbNvG5/PR2tra7+Mx02aNx+JPTRaWDSUuuH66k7MnGLrdjQwb\nBSwREclapmnS3NxMZ2fnMY/Zts2WDpvfHjLpjIMBnF/t4APTHBoaKsNOAUtERLJSPB6nqamJQCBw\nzGMHAha/O2RxMJg6HTi5GD5S52SGW3OtZGQoYIlkmYRpkTBtkpZNwrJJmBZJy8aywDDAYQAYOIzU\n93lOBwVOgwKXA5dDl5/L2BCNRvF4PEQikT7HfRGb3zeYbD0y06rUBVdPcXBejQOnfvdlBClgiWRY\nIJakJRCnI5KkM5KgI5ygI5ykI5IgEEsSipuE4hahuEk4bpKwBr6P2kAcBuQ7HRS6HLgLnJQWOCnN\nd+EucB75cFF65HNFkYvxxXmML87DXeBUMJNRIxQK0djYSCKR6D3WFUv1Wf3VZ2GRujrw4loHl05y\nUOjU766MPAUskREQSZgc6oxysDOKxx+jyR+jKRCjyR8nGO//cvLj6QlJeU6DPIeBy2ngcjhwGmAD\nlp3qPUl9bZMwbWJJi1jSwrQhmrSIJi26osefEfReeQ6DyuI8xpfkMb4oj8riPKpK8qhx5zPRnc9E\ndwFlCmEyAvx+Px6Pp3cMQ3vUZp3X4s3W1O+3AZxXbXD1FCfl+fp9lMxRwBJJs45wgh2+EHvawhzs\njHKgI0JTIH7c9YUuBzXufMYX51FZ5KKyOI+Kojwqi12UFbooyXNSkp/6KM53UuA88Wk+27aJRqNE\nIhGi0SjBYJB4PI5p2SQsiFsQTkLEtAkne76GqGkQsx1ELYNAArqiFp1Rk1DCoiUYpyV4/D9HUZ6D\nmtJ3A9dEdz6TywqYUl7IRHc+Tof+sZOh6ezsxOv1pq4ajNis9ZpsaLWxSAWrReMNrprsZGKxftck\n8xSwRIbAtGz2tofZ4QuxvSX1ub8Q4nIYTBtXwIyKIqaUF1BbVkCtu4BJZfmUF7qGvPNjWRaxWIxI\nJEIkEiEYDPY5fdLD6TBwOqAQKMuH1D9L72Uf+eh9FgnLSchyErZdhEwngSR0xmzaIhatoSQtoQSh\nuMnBI7t0/f35J5UVMKW8gKnlBUwZV3jk60LKCvXXkAzMtm1aW1tpaWnhYNBmfbPF5vbULq0DOHuC\nweWTndQUKVjJ6KG/2UROgW3beP0xNngC1HsCbG4KEnrPKb7iPAdzq0qYU1XMzMpC6iqLmFJeiCuN\nOziWZfXuUEUiEUKhUL+BKl3yHDDOYTIOE/KOHCxNfTIMA5ergKQjD7/poiNm0xGD1rBJUzCJJxCn\nLZSgoStKQ9ex4auswMmUcYVMH1fI9IpCZlQUMr2iiMqioQdPyX6madLobWbt3nZebTZpDKWOOw04\np8rgsklOJmhQqIxCClgiJxA3LTZ5g/zlUBcbGgPH7FBNKitgfk0J82pKOKO6hGnjCtN+OsyyrN7T\nfSMRqE6FbdtHakngBtxOmF4MFIOj2oHTmYeRV0pH3IEvmjq10xxM4g0m8Pjj+GMm21tCbG8J9Xld\nd4GTGRVFRwJXKnjNqCjSjlcOaewM8Wz9YV4+HCF0pGWwxAXnVTu4sMbBON03UEYx/U0l0o9owuSt\nxgCvHezirw3dhBNW72NlBU4WT3KzZLKbJZPLqHHnp/39TdM8ZodqoBvXjlaWZWFZFiQSlJK6ZH6m\nG3CDw+HA4SggZhTQFjdojtg0BU0aA0kaumMEYiZbm4NsbQ72ec3KIhfTjwSvGRWFzKgsYvq4Qorz\nnRn5M0p6RRIm6w90sXpXG1tbwr3HJxfDRbVOFo9PXdwhMtopYIkcYVo2bzf6eWlPB39t6CZuvtuH\nNLOykAtnjOPcaeXMHl+U9ttrmKZ5zA5VNgaqU9ETvpwkqTGgphgWFgPVYBguQpaT9oSL5gh4wxae\nQJJGf2qcRUckwEZv3+GSNaX5fULXjIoipo4rIN+pwZKjXdK02OTxs2ZvB68d8hM78r89lwFnVRos\nm+hgRqlmuEl2UcCSnLe/I8JLuztYt6+Dzsi7oeaM6hIumFHOhdPHMbm8IK3vmUwmiUajRKNRwuEw\noVCo97JzSZ12LDaSFOcnmZoPlAO1YNlOuhMOWuMOfDEHzWELT9DEE0j0XuX4xmF/7+s4DJhcnrq4\noO7IKcYZlYXUugt0VWOaWZaFaZq9wbnno+fY0Y+Zpkk0nuSd1ij1vgRb25O9pwAB6twG51Q5WFhp\nUKRb2kiWUsCSnBRNWvxpbwe/29HG3vZ3J0FPLS/gytMquWx2JdWl6Tv11xOoIpEI4XCYcDisQDUI\nDsOgIt+mIt9kTqkJ41PHTdtFWxR8UYOWqEFzFLxBk5aQyeGuGIe7Yqw/8O7r5DsNpo1LXYAwvSd4\nVRRSVZKX07sk/QWhgYKTaZokk0mSySSWZaXmrx356I8/brO722Z7l8WOTpvYu2feqSqEReMdnFPl\nUNO6jAkKWJJTWgJxfrejlT/saicQSwWc0nwnl86q4MrTKplbVZyWf2ATicQxO1SWZZ34iTIoTsOg\npghqiuCs3hETDhKWQUsEmsI2LVF6v+6IWuxtj/QJ1wAl+c7e04zTK4pSIyXKC6kqzUv7aeHh0hN8\njg5Ctm0fNzgdHZROJiSdiphpcyBgs6vbZne3RVO47+OTi2FBpYOzKh3UFJHT4VbGHgUsyQlbmoL8\napuP1w9103OnmblVxXzozCoumjGOfNfQ+nR6AtXRO1QKVJmX5zCYUgJTSvr+wx1JOlJN9WGb5gg0\nh6EpYhOKm2xrCbHtPVc0FjgNJh8JW1PKC5h6ZI7XlPJCStLcXP/eMNTfrtLRa3rC0dHBKp0h6VTq\n7orDwYDNoaDNwaCNJ2RzVCsjeQ6Y5TaYO87grEoHlboKUMYwBSwZs2zbpt4T4MlNzWxtTv2D6XIY\nXDJzHB86s4rTq0sG/drxeLxPoIpEIgpUWaTIZVDnNqhzv3vMtm0CCXqDV0vEpjWaunlwIGGzvyPK\n/o5j53hVFLmYUl5I7ZHbBtW4C6gpzaO62MW4Qgf2e06t9Xfq7b07ST3haLT+Tpm2TVs0tRvY89EQ\ntPG/Z3KIAUwtMZhTnvqocxtpnQcnMpopYMmYY9s2bzT4eWJTM7taU+ck3AVObjijiuvnTaCyOO8E\nr3CsowNVKBQiGo2O2n/8ZHAMw6AsH8ryDeaU930skrRpjdr4IuCL2rQeCV+tUZvOSJLOSJCtzce+\nptOAygKoKDCoOPLa5flQlpd6r/J8g9K81CnO0caybbrj0B6zaY+m7vnXHkv9d2gOQ7KfzbEiJ8xw\nG0wvTX1MK1WTuuQuBSwZM3qC1U83NLG/I9VbU17oYsVZ1Vw3b8IpncqJx+N9ZlDFYjEFqhxW5EqF\nhWmlfY9bR06LtUZS0+s7YqkQ0hlLBZNggt4glnJsKjGA0jwoy4Nil0GRC4pdR752pj4Xu6DIBfmO\n1GnPPEfqdFu+A1yO1NWSPa9lGKnPtp0KQaYNSevdz1EToqZN1Dxy/8mkTSgJgYSNPw6BxJGvE/Se\nTu9PRT7UFhu9H5NLDKoKyZpeNZHhpoAlY8JOX4gfventHUpZWezixgU1XDt3PIV5Jw5WPffxi0aj\nvTtUI9m/ItnJYRhUFnDcXqKY+W7Y6o6nAoz/SJDpjqdCTDDRE2qgbwDL/O+fOw/GFxiML+z5bDCh\nECYWaWdK5EQUsCSref0xfvKWl5cPdAGpKesfWzyRD5w+YcDG9aMDVTAYJBaLKVBJ2hU4DSYWw8Ti\n44cR07J7d43CydSuUjh55Otk6uvUMYhbkLBsEhYkLIibqZ0py373Ft32ka8dgNORGtbZ+9mAQmdq\nR67Q+e7XRc7Uzb/deQbuvNSpTHcempguMgSDDlj/+7//y969e3E4HNTU1HDHHXfw+uuvs3r1ahwO\nB3V1ddx+++3prFWkVzCWZFV9M7/b0UbSssl3Gnx4fjU3Law55lSgbdvEYrE+PVQKVDJaOB0G4wrQ\nffVExphBBaxwOMzu3bv5+te/DsC3vvUtfD4fa9euZeXKlRiGwbPPPsu6deu47LLL0lqw5Dbbtlm3\nr5Mf/tVDVzSJAVw9p5JPLK2lqiS/d03PDKqeHap4PK5AJSIiI2ZQAau4uJgLLriAz372s+Tl5bF0\n6VK2bdvG8uXLewfFXXHFFfzoRz9SwJK0aeiM8r2/HGZzU6rP6syaEu66YCozKwuJRqN0dnYSiUQI\nBoMkEgkFKhERyZhBBay2tjbefvttHn74YfLz8/n+979Pd3c3M2bM6F1TVlaG3+8//ouInKRowuSJ\nTS08u6UF04byQiefXFTNebV5xGOd7NmT2qESEREZLQYVsPbs2cPixYspLCwE4MILL+TRRx/llltu\n6V3j9/spKytLT5WSs7Y2B/n2y4doCsQxgEumFnLNZIN82mhuynR1IiIi/RvU/UGmTJnCtm3beucC\nbdmyheuvv57169f3HluzZg1Lly5NX6WSU2JJi0f/2sg/Pr+HpkCcScUGXzzTyQ2TTfJJZro8ERGR\nAQ1qB2sSqvG5AAAgAElEQVTq1KksXbqU+++/H8MwmDdvHjfccAPjxo3jvvvuw+l0Mn36dFasWJHu\neiUH7PSF+NbLhzjcHcNhwFVTXFxei26xISIiWWPQYxquuuoqrrrqqj7Hli9fzvLly4dclOSmhGmx\nqr6ZX2xpwbJhankBHzvNRbVL/VUiIpJdNGhURoXmQIyvrzvIrtYwBvDRs6q5stYm6O/KdGkiIiKn\nTAFLMm79gS6+u76BUNykpjSfL18ynRpXlObmfu6eKyIikgUUsCRj4kmL/37Tw2+3twFwwfRyvrR8\nGsQjNDS0ZLg6ERGRwVPAkozwdMf4+roD7G2P4HIY3HHuJD54RhXRaJRDHo+GhIqISFZTwJIR9+bh\nbr6x7iDhhEWtO59/vqyOOVXFJJNJmpqaSCY1hkFERLKbApaMGNu2+cUWHz9+y4sNLJtRzpeWT6ck\n34lt2zQ3NxMOhzNdpoiIyJApYMmIiCUt/n19A+v2dQLwyaW13LKoBseRe1e2tbXR1aUrBkVEZGxQ\nwJJh1xqK87WX9rOnLUKhy8E9l0znwhnjeh/3+/34fL4MVigiIpJeClgyrHb6Qnz1pf10RpJMdOfz\nL1fOpK6yqPfxSCSCR03tIiIyxihgybD566Fuvr7uADHTZmFtKfdeXkd54bu/colEAq/Xi2maGaxS\nREQk/RSwZFj8fmcb33vtMJYNV8+p5P8sm9bnXoKWZdHc3EwkEslglSIiIsNDAUvSyrZtVtU38/jG\n1BT2jy+eyG1LJmIYfW/U3NbWRnd3dyZKFBERGXYKWJI2pmXz8GuHeXFXOw4DvnjhVD5w+oRj1nV3\nd6upXURExjQFLEmLaNLi62sP8MZhPwVOg69cVsf508uPWRcOh/F6vRmoUEREZOQoYMmQRRIm9/9x\nP5ubgpQVOFl59SzmVZccsy4ej+PxeNTULiIiY54ClgxJKG7yz6v3sb0lRGWxi2++/zSmVRQes66n\nqT0Wi2WgShERkZGlgCWD5o8m+efV+9jVGqaqJI9vXnsak8sL+l3b2tqK3+8f4QpFREQyQwFLBqUr\nkuD/vbiPfe0Rat35PHTtbCa6+w9XnZ2dtLa2jnCFIiIimaOAJaesI5zgnhf2cqgrypTyAh66djZV\nJfn9rg2FQjQ1NY1whSIiIpmlgCWnpCuS4J9e2MPhrhjTxxXy0LWzqSzO63dtPB7H6/ViWdYIVyki\nIpJZClhy0vzRJF/+wz4Od8Woq0iFq3FF/Ycr0zRpampSU7uIiOQkR6YLkOwQPnK14P6OCFPKC/i3\nAcIVgM/nIxAIjGCFIiIio4cClpxQNGlx3x/3s6s1TE1pqqG9YoBw1dnZSXt7+whWKCIiMrooYMmA\n4qbFA2v2s7U5yPjiPL45QEM7QDAY1KR2ERHJeQpYclxJy+Yb6w7ydmOA8kIXD107m9qy/kcxAMRi\nMTweD7Ztj2CVIiIio48ClvTLtm0efrWBvxzqxl3g5N/eP5tp446d0N7DNE28Xi+JRGIEqxQRERmd\nFLCkX6vqm1m9u4MCp8HXr57FrPFFx11r2zYtLS2EQqERrFBERGT0UsCSY/xhZxuPb2zGYcA/X17H\n6f3cuPloHR0ddHR0jFB1IiIio58ClvTx5uFuHn7tMABfvGAq500rH3B9IBCgubl5JEoTERHJGgpY\n0mtXa4iVaw9i2XDroho+MG/CgOuj0aia2kVERPqhgCUAeP0x7l29n1jS4srTKvnk0toB1yeTSbxe\nL8lkcoQqFBERyR4KWII/muSfX9xHdzTJkslu/u+yqRiGcdz1PU3t4XB4BKsUERHJHgpYOc48MuvK\n448xs7KQ+y6vI8858K9Fe3s7nZ2dI1ShiIhI9lHAynE/fMNDvTc1SPRfrpxFSb5zwPV+v5+WlpYR\nqk5ERCQ7KWDlsD/sbOPX21pxOQy+dkUdNe7j3wIHIBKJ4PV61dQuIiJyAgpYOWprc5Dv/aURgL+/\ncCpnTiwdcL2a2kVERE6eAlYOagnEeWDNAZKWzYfnV3HN3PEDrrdtm+bmZiKRyAhVKCIikt0UsHJM\nJGHy1ZfevWLwM++bfMLntLW10dXVNQLViYiIjA0KWDnEtm2+80oD+zuiTCkv4J8vm4HTcfxxDJBq\navf5fCNUoYiIyNiggJVDfr2tlVcOdFGc5+BrV87EXeAacH0kEtGkdhERkUFQwMoR21tC/PcbHgD+\nYfk0po0rHHB9IpHA4/FgmuZIlCciIjKmDLyFcQI+n49Vq1Zx991343A4WL9+PatXr8bhcFBXV8ft\nt9+erjplCLoiCb6+7gCmDR+eX8XyuooB11uWRXNzM9FodIQqFBERGVsGvYNlWRbPP/88d955Jw6H\nA5/Px9q1a1m5ciUPPPAAbrebdevWpbNWGQTTsvm3Px+iNZTgjJoSPn2STe3d3d0jUJ2IiMjYNOiA\n9Zvf/IZDhw7x8MMP86c//YnNmzezfPny3nvYXXHFFWzYsCFthcrgPLGxmXpPalL7vZfNwHWCpvau\nri41tYuIiAzRoE4R+nw+Ghoa+OpXvwrA97//faZMmcK0adN615SVleH3+9NTpQzKW4f9PLGxGQP4\nf5dOZ0LJwJPaw+EwTU1NI1OciIjIGDaoHaz6+nrOOeccHA4HDoeDZcuWYVlWn0Dl9/spKytLW6Fy\nanzBOA/9+SA28ImltSyZPPDPIh6Pq6ldREQkTQYVsNxuN1u2bOn9ftOmTRQWFrJ+/XosywJgzZo1\nLF26ND1VyikxLZtvvnwIf8zknCll3LKoZsD1PU3tsVhshCoUEREZ2wZ1ivCCCy5g37593HvvvRiG\nwfz587nuuusoKyvjvvvuw+l0Mn36dFasWJHueuUk/GJLC1uaglQWufini6fhMAbuu2ptbdXpXBER\nkTQaVMAyDINPfOITxxxfvnw5y5cvH3JRMng7fSH+d0Oqj+qfLp7OuKK8Add3dnbS2to6EqWJiIjk\nDA0aHUPCcZMH/3QQy4aPzK9m6ZSB+65CoZCa2kVERIaBAtYY8sjrjTQF4swaX8Tt59QOuDYej+P1\nent75kRERCR9FLDGiHV7O3hpTwcFToP/d+kM8p3H/9GapklTU5Oa2kVERIaJAtYY0ByI8Z+vHQbg\n8+dPOeF9Bn0+H4FAYCRKExERyUkKWFnOtGwe+vMhwgmLC2eU8/654wdc39nZSXt7+whVJyIikpsU\nsLLcr7b52NYSYnxxHncvm9Z7q6L+hEIhvF7vCFYnIiKSmxSwslhDZ5SfvJ26CvDui6ZSVnj8qRux\nWIzGxkZs2x6p8kRERHKWAlaWMi2bb71yiIRpc/WcSt43tfy4a5PJJF6vl0QiMYIVioiI5C4FrCz1\nzFYfu1rDVJXk8bnzphx3nW3b+Hw+QqHQCFYnIiKS2xSwstDBjgirjkxr/4eLplGS7zzu2o6ODjo6\nOkaqNBEREUEBK+skLZtvvXyIhGXzgdPHDzitPRAI0NzcPILViYiICChgZZ2nN7ewpz1CTWk+n3nf\n5OOui0ajeDweNbWLiIhkgAJWFtnXHuHx+tSpwS8tn0bxcU4N9jS1J5PJkSxPREREjlDAyhKmZfPd\n9Ycwbbh+3gQWTXL3u862bVpaWgiHwyNcoYiIiPRQwMoSv9neyp62CFUleXzqnEnHXdfe3k5nZ+cI\nViYiIiLvpYCVBXzBOD89MlD0ixdMPe6pQb/fT0tLy0iWJiIiIv1QwBrlbNvme68dJpq0uKhuHOdN\n73+gaCQSwev1qqldRERkFFDAGuXWH+jijcN+SvKdfOH8/geKqqldRERkdFHAGsUCsSSPvN4IwKfP\nmcT44rxj1ti2TXNzM5FIZKTLExERkeNQwBrF/uctL52RJGfWlPD+08f3u6atrY2urq4RrkxEREQG\nooA1Sm1tDvLCznZcDoP/u2wqDsM4Zk13dzc+ny8D1YmIiMhAFLBGoYRp8R+vNgBw88IaplcUHbNG\nTe0iIiKjlwLWKPSrd1o53BVjSnkBNy+sOebxRCKBx+PBNM0MVCciIiInooA1yrSG4jy+MXWD5jvP\nn0K+q++PyLIsmpubiUajmShPREREToIC1ijzoze8RJMWy2aUs3RK2TGPt7W10d3dnYHKRERE5GQp\nYI0im7wB/ry/kwKnwWfPPXbmVVdXl5raRUREsoAC1iiRtOzemVc3L5pIjTu/z+PhcJimpqZMlCYi\nIiKnSAFrlPjNtlYOdUaZVJbPR8+q7vNYPB5XU7uIiEgWUcAaBdrDCVbVp3anPn9e38b2nqb2WCyW\nqfJERETkFClgjQL/86aHcMLi3GllnDut782cfT4ffr8/Q5WJiIjIYChgZdg7zUHW7O0kz2nw+fP6\nNrZ3dnbS1taWocpERERksBSwMsiybf7rSGP7jQtqmFRW0PtYKBRSU7uIiEiWUsDKoDV7OtjbHmFC\ncR43LXi3sb2nqd2yrAxWJyIiIoOlgJUh0YTJT99O7VDdfs4kCvOcAJimidfrJR6PZ7I8ERERGQIF\nrAx5dquPtnCC08YXcfnsit7jPp+PYDCYwcpERERkqBSwMqA9lODpLamJ7J89bzIOwwCgo6OD9vb2\nTJYmIiIiaaCAlQE/3eAllrS4cEY5C2rdgJraRURExhIFrBG2rz3MH3d34DTg0+dMAiAWi9HY2Iht\n2xmuTmT4/fKXv+Saa65h+/btmS5FRGTYKGCNINu2+eEbHmzgg2dWMbm8kGQyidfrJZFIZLo8kRFx\nzTXXUFBQwLx58zJdiojIsFHAGkFvNPjZ5A3iLnBy66KJ2LaNz+cjFAplujSREVNfX8/ChQsxjvQe\nioiMRQpYIyRp2fzoTQ8AH188kbJCFx0dHXR0dGS4MpGR9dZbb2EYBi+99BLf/va32bdvX6ZLEhFJ\nOwWsEfLirnYOd8eYXFbAdfMmEAgEaG5uznRZIsPq2Wef5brrruPTn/40hw8fBmDDhg3ceOONXHnl\nlVx44YU89thjGa5SRCT9hhSwYrEYX/3qV1m1ahUA69ev59577+X+++/nJz/5SVoKHAuiSYvH63uG\nitZiJlKT2tXULmNZfX09jzzyCA899BDhcJiHHnoIn8+HZVmceeaZALS3t9PV1ZXhSkVE0m9IAevx\nxx/n0ksvBaC1tZW1a9eycuVKHnjgAdxuN+vWrUtLkdnu19ta6YgkmTOhmPOnlOL1ekkmk5kuS2RY\nPfroo7zvfe9j9uzZ2LZNdXU1O3fuZOHChb1r3nrrLc4777wMVikiMjwGHbD+8Ic/sHjxYqqrU/fQ\n27RpE8uXL+9tXL3iiivYsGFDeqrMYoFYkl9sbgHg9rNr8fl8hMPhDFclMry2b9/Orl27uOyyyygo\nKOCpp57i/vvvp6SkhNLSUgAOHz7M/v37uemmmzJcrYhI+g0qYO3YsQO/38+SJUt6T3MFAgHKysp6\n15SVleH3+9NTZRb7xeYWgnGTxZPcTC+M09nZmemSRIbd6tWrATj33HP7HF+yZAmGYfDiiy/y7LPP\n8h//8R8UFhZmokQRkWHlGsyTNm3aRENDA9/61rcIBAJ0d3dz0UUX9QlUfr+/T+DKRW2hOL/a1grA\nTWeU0dLSkuGKREbGa6+9Rl1dHePGjetz3DAMvvCFLwCpeVgiImPVoALWLbfc0vv19u3b2bBhAxdf\nfDH/9V//xSWXXILD4WDNmjUsXbo0bYVmoyc2NhM3bS6Y6qY42kFSTe2SAw4fPkxrayvnn39+pksR\nEcmYQQWs9zIMg6qqKi699FLuu+8+nE4n06dPZ8WKFel4+azk6Y7yh13tOAy4stZSU7vkjPr6egDO\nOOOMDFciIpI5Qw5YZ5xxRu9fpMuXL2f58uVDLmos+OmGJiwbLppSSLkjnulyREbMxo0bAZgzZ06G\nKxERyRwNGh0Ge9vCvLy/izwHXFKlnSvJLZs2bSI/P58ZM2ZkuhQRkYxRwBoGP3nbC8CFNQ4qCnS/\nNckdhw8fprOzk7q6OpxOZ6bLERHJGAWsNNvWEuStxgAFTrhskv7zSm7ZvHkzAKeddlqGKxERySwl\ngDRbVZ+6v+BlU/IpzdPuleSWnoA1a9asDFciIpJZClhptK05SL0nQHGeg8unFmS6HJERt23bNgBm\nzpyZ4UrANM1BP1dX/YrIUClgpVHP7tWHzqyiRLtXkmO6urrweDwYhkFdXV1Ga3n55Zf54x//OOjn\nP/7447zzzjtprEhEco0CVpq80xyk3pvavfrIWdWZLkdkxPUEknHjxlFeXj7s79fY2Mg999zDo48+\nyne+853e23Zt2rSJzZs38/73v3/Qr/3xj3+cn/3sZxw8ePCkn/PII4+wYsUKLr74YjZt2jTo9xaR\nsUEBK016dq8+PL8ad0Fa5reKZJWegDUS/VeJRIJ//Md/5JJLLqG9vZ3f//73hEIhQqEQP/jBD/js\nZz87pNd3uVx86Utf4hvf+MZJny688847ufXWW8nLy9OQVRFJzyT3XLe1OchGb4CSfCd/M78q0+WI\nZMT27duBkem/evPNN2lqamLhwoXMmDGDa665htLSUh599FGuvPJKCgqG3gNZU1NDXV0dL774Itdd\nd91JPWfLli2cfvrp5OfnD/n9RSS7aQcrDVbVNwHw4flV2r2SnGSaJrt27QJGZgdr06ZNjBs3jkmT\nJjFv3jyWLl1KJBLh+eef5+qrr07b+3zkIx/hySefPOn1W7ZsYdGiRWl7fxHJXgpYQ7SlKcgmbzC1\ne3Wmdq8kNzU0NBCNRjEMg9mzZw/7++3YsYPTTz+9z7HXX3+d2tpa3G532t5n9uzZdHd3s3v37hOu\n9Xg8tLe3K2CJCKBThEN29O5VqXavJEft3LkTAKfTOay3yPnGN75BZ2cnW7duZdq0afzTP/0TkyZN\n4u677+btt99m/vz5x33url27WL16NU6nk6amJu655x5++9vfEggEaGtr4+/+7u+YNGlSn+c4HA4W\nLFjAW2+9dcy9FTds2MDvfvc7Jk6cSDAYZObMmTidTs4666whv6+IZD8lgiHY0hRgc5N2r0R6AtaM\nGTNwuYbvr5WvfOUreL1ebrnlFu644w4uuuii3sf27t3L9ddf3+/zvF4vL7zwAnfffTeQCmqf+9zn\n+MpXvoJt29x1113MmTOHG2+88ZjnTpkyhb179/Y59vzzz/PYY4/x2GOPMWHCBHw+Hx/72Mc4/fTT\n+/R/DeV9RSS76RThEDy5sQXQ7pVIT8B67y7PcNizZw/AMacim5qaKC0t7fc5Tz/9NJ/73Od6v49G\no5SVlXHmmWdSXV3NTTfddNyxDm63m6ampt7v9+7dy3e/+13+/u//ngkTJgBQXV1NYWHhMacHh/K+\nIpLdFLAGaacv1Dv36kPavZIcZpom+/fvB0YmYO3du5eSkhJqa2v7HA+FQscNWDfffDNFRUW932/b\nto2lS5cCqXD0+c9//ri9W+Xl5YRCod7vf/SjH1FSUsLFF1/ce+zgwYP4/f5jAtZQ3ldEspsC1iD9\nfHNq9+r6eRN05aDktIaGBuLxOIZhHNN4Phz27t3b782kDcPoHTb6XkeHsYaGBtra2li8ePFJvZ9l\nWb2vGwgEePPNNzn77LNxOp29azZu3IjD4Tim/2oo7ysi2U0BaxAOdkb4y6Fu8p0GH9bUdslxPf1J\nLpdrRK4g3Lt3b7/vU1pait/vP+Hz6+vrycvL69MQ7/V6j7s+EAj07ox5PB4sy+LMM8/ss2bjxo3M\nnTuXoqKi477Wqb6viGQ3BaxBePrI7tXVc8ZTUZSX4WpEMqsnYM2cOZO8vOH934Pf78fn8/UbsGpr\na/sNWLFYjB/84Ae9pzHffvttZs2a1duMblkWTz311IDv2XOVX0lJCZAaQnr062/atKn39OCzzz6b\nlvcVkeymc1unqMkf40/7OnEY8NEF2r0S6QkQ8+bNG/b36mlw72+Y6VlnndXvvQP/+te/8vOf/5y5\nc+fidDppbGzs06u1atWqARvNDx48yNlnnw2kriicNWtWb9N7Mpnku9/9LslkkkmTJtHZ2UllZWVa\n3ldEspsC1il6ZqsPy4YrT6tkonvot+MQyXYjGbB2795NaWlpvztY5557Lt/73veOOb5o0SLe//73\ns2vXLnbv3s2jjz7Kd7/7Xb797W+Tl5fHsmXLjnvvwGQyyTvvvMPnP/95INXn9S//8i98//vfx+fz\nYVkWn/jEJ1i8eDEvvvgiu3fv7r0P4lDeV0SynwLWKWgPJ1i9qx0DuGlBzQnXi4x1PcMyDcMYkbCw\ne/duli5disNxbHfDggUL6OjooK2trXd8AqSuAvzyl7/cZ+1XvvKVk3q/nTt3UlNT02fHbOrUqTz0\n0EN91k2ePJlrrrmmz7GhvK+IZD/1YJ2CX271kbBsLpxRzrSKwkyXI5JxPafk3G4306ZNG5b3eOKJ\nJ/jSl74EpALPJZdc0u+6/Px8PvzhD/f2QKXDL37xC2666aa0vZ6I5A4FrJPkjyZ5fmcbADcvnJjh\nakRGhwMHDgCp3aPh8tJLL5GXl8e+fftwuVx95k+91y233MJf//pXAoHAkN+3oaGBlpYW9UmJyKAo\nYJ2k325vJZKwWDrZzZyq4kyXIzIq9OxgLVy4cNje4+abb2bChAmsWrWKf/3Xf+0zf+q9CgsLueee\ne/jmN7953JlYJyMWi/Hwww9z//33YxjGoF9HRHKXerBOQiRh8qttrQDcski9VyI9ekY0vHfAZjpd\nc801x/Q3DWTevHnccMMNPPfcc6xYsWJQ7/n444/zmc98hsmTJw/q+SIiClgnYfXuDgIxk3nVxZw1\nsf9bcYjkon379lFYWDgit8g5Feeccw7nnHPOoJ//qU99Ko3ViEgu0inCEzAtm1++4wPgo2fV6HSB\nyBHNzc0Eg0HOOOOMAU/biYjkIgWsE3jtUBfNgTiTyvI5f3p5pssRGTV27doFoHvriYj0QwFrALZt\n88yW1O7Vh+dX43Ro90qkx44dOwBYunRphisRERl9FLAGsK0lxK7WMGUFTq6aMz7T5YiMKtu3b6ek\npGREJriLiGQbBawBPLs1tXt13bwJFLr0n0qkRywWY/v27Zx99tn9TlUXEcl1+pvxOBq7o7x+qJs8\nh8EHz6zKdDkio0p9fT2JRIJly5ZluhQRkVFJAes4fvlOKzZwxWmVVBTlZbockYz6z//8T26//XaS\nySSQmq7udrsHnKouIpLLFLD60RVJ8Mfd7QB8ZH51hqsRyby3336bWCyGZVn4fD5eeeUVbrzxRgoK\nCjJdmojIqKRBo/14fkcbcdPm3KlluqmzCKlJ7RUVFQQCAR588EGmTJnCrbfemumyRERGLe1gvUcs\nafGb7ambOq9YoN0rEYA77riDHTt2cMstt5Cfn8+3v/1tXK7+//9ZMpnkscce4ze/+Q3PPfccX/7y\nl/F6vSNcsYhIZmkH6z3W7O2gO5rktAlFLNBtcUQAKC8v5zvf+c5Jrf3Od77DrFmz+OAHP0h3dzc/\n/vGPmTRp0jBXKCIyumgH6yi2bfOrd1I3dV5xVrVuiyNyivbt28e6deu44YYber9ftGhRhqsSERl5\nClhHqfcGaOiKMqE4j4vqKjJdjkjWefvtt1m4cCH5+fkAbNiwgSVLlhAIBDJcmYjIyFLAOsqvt6V2\nr64/YwIu3RZH5JSVlZVRWVkJQDgc5pVXXmH+/Pm89NJLGa5MRGRkKWAd4emO8WaDnzynwfvn6rY4\nIoNx+eWXYxgGa9eu5eWXX+bKK6/kjTfeoLa2NtOliYiMqEE3ub/wwgu8+uqruFwuamtr+cxnPsPr\nr7/O6tWrcTgc1NXVcfvtt6ez1mH12+2pwaKXzapgnAaLigxKfn4+99xzT6bLEBHJuEEFrGAwSEND\nA1//+tcxDIPHH3+c3//+92zcuJGVK1diGAbPPvss69at47LLLkt3zWkXjpusPjJY9EO6LY6IiIgM\n0aBOEZaWlvK5z32u9yq7WCyGbdssX76899gVV1zBhg0b0lfpMHppTwfhhMVZE0uZNb440+WIiIhI\nlhtyD9avfvUrioqKsCyLsrKy3uNlZWX4/f6hvvyws2yb32xPNbdr90pERETSYdABy7IsHnvsMfLy\n8rj11ltxu919ApXf7+8TuEarDY0BGrtjVJXkccH08kyXIyIiImPAoAJWNBrl3//931mwYAHXXXcd\nAAsXLmT9+vVYlgXAmjVrWLp0afoqHSY9oxluOKMKp0YziIiISBoMqsl93bp17N27l2AwyB/+8AcA\nLr30Ui699FLuu+8+nE4n06dPZ8WKFWktNt0au6O81egnX6MZREREJI0GFbCuvfZarr322n4fW758\n+ZAKGkm/2Za6qfPlsyspK9RtGUVERCQ9cnbQaChu8sc9qdEMH1Rzu4iIiKRRzgasl/Z0EElYLKgt\nZWZlUabLERERkTEkJwOWbds8vyN1evCGMyZkuBoREREZa3IyYG1tDtHQFaWyyMUF08dluhwREREZ\nY3IyYPXsXl0zdzwujWYQERGRNMu5gNUZSfDqwS4cBlx7uk4PioiISPrlXMBavaudpGVz7tRyqkvz\nM12OiIiIjEE5FbBMy+b3O1OjGa6bp90rERERGR45FbA2ePy0BONMdOezdIo70+WIiIjIGJVTAet3\nR5rbP3D6BByGmttFRERkeORMwGoJxHmzwU+ew+DqOZWZLkdERETGsJwJWC/sasMGLqobx7iivEyX\nIyIiImNYTgSshGnx4q5Uc/sH1NwuIiIiwywnAtZfDnXTGUkyvaKQ+TUlmS5HRERExricCFg9k9uv\nnzcBQ83tIiIiMszGfMBq7I6yuSlIgcvB5bPV3C4iIiLDb8wHrNVHeq8umTmOknxnhqsRERH5/+3d\nzU8UeR7H8U836CLSzcoMoK5AnAVHDRpp0CVoyLLozGxiTDQe9MDBgw+JXjwQL6jRjrqJB4N/gInR\n6NFgjCwz2m6Cm834wIMsEAHHofch2NuLQ/MgKNO9B4aeNZmgjWVVV/X7dasEqr7wbeBD/aq/P6QC\nR++sOk0AAAahSURBVAes6WhMX/cPS5rZ2BkAAMAMjg5Y3wZnHm4v/HWG1ubxcDsAADCHowPW7GiG\nP37+CQ+3AwAA0zg2YIXHX+vhPyNKd7tUW7zE6nIAAEAKcWzA+rpvWNGYVFWUzeR2AABgKkcGrGgs\npj/3/bw8CAAAYCZHBqzOf49paPS18rMWquw3HqvLAQAAKcaRAav5p4fbv1iVIzcPtwMAAJM5LmBF\nJqf11+9/kEvSl6tYHgQAAOZzXMC6OzCsN9GYKlZ4lJe10OpyAABACnJUwIrFYvHlwa8+/9TiagAA\nQKpyVMB6+p8Jff9yUtkZ6aos9FpdDgAASFGOClizoxm2leRoQZqjvjQAAGAjjkkhU9NR/eXZS0nS\nl6tyLK4GAACkMscErL8NjmjiTVSrPs1U0ZJFVpcDAABSmGMC1jf9P8++AgAAsJIjAtZ/x9/o8b9G\nle526fefsbEzAACwliMC1t2BmY2dKwu98makW10OAABIcbYPWLFYTN/0D0uSvihhcjsAALCe7QNW\nX3hCgz/MzL6qKGD2FQAAsJ7tA9bs3ava4iVKd7OxMwAAsJ6tA9brH6O699Psq20sDwIAgCRh64D1\nbTCi0akf9VnOIv32E2ZfmWloaMjqEjBP9M7e6J990bvUYuuANTv7alsJs6/M9uLFC6tLwDzRO3uj\nf/ZF71KL4TMNbt68qYcPH0qSfD6fdu7cafQlJEkvX73Rg39E5HZJfyhm9hUAAEgeht7B6u3t1fPn\nz+X3++X3+zU0NKSuri4jLxEXGHipaEzaVODVkkULPso1AAAA5sPQgNXe3q7a2tr4cW1trdra2oy8\nRNzsuwd5uB0AACQbQ5cIR0dH5fF44sder1cjIyNGXkKS9N3wK303/EqeX6Xpd4XJOftq8eLFWrDA\nuXfWcnJy5PUm5/cec6N39kb/7Muq3i1cuND0a8LggOXxeBSJROLHkUhkzhfTh9zd+pNPkqb1986O\neZ8D85efn69wOGx1GZgHemdv9M++rOpdOBzW4OCg6ddNdYYGLJ/Pp5aWFq1bt06SFAgEtHnz5l/8\n2P9fSgQAAHASVywWixl5wqamprfeRbhr1y4jTw8AAJD0DA9YAAAAqc7Wg0YBAACSEQELAADAYAQs\nAAAAgxGwAAAADGb4XoTvYtZehfg4bt++rfv37ys9PV3Lli3T/v37lZ5u+ssIH2Bqakpnz55VcXGx\n6urqrC4HCQiFQrpy5YqOHj0qt5v/j+3k8uXLGhgYkNvtVn5+vg4cOMDvziQ2PT2ta9euqbe3V+fO\nnZMktba2qqWlRW63WytXrtS+ffvmPIepP6Fm7lUI442NjSkYDOrMmTM6ffq0PB6PHjx4YHVZSNDV\nq1dVU1NjdRlIUDQa1a1bt3T48GHClc1MTEyor69Pfr9fp06d0vj4uEKhkNVlYQ7Xr19XaWlp/DgU\nCunu3bvy+/3xv3+BQGDOc5j6U2rmXoUwXlZWlg4dOiSXyyVp5k5Ibm6uxVUhEc3NzSorK1NeXp7V\npSBBTU1NGhwcVGNjo+7du2d1OUhAZmamqqqqdPDgQR05ckS5ublavny51WVhDnV1dfL5fPHjjo4O\nVVdXx//+bd26VY8fP57zHKYGLLP2KsTHd+PGDWVmZqqkpMTqUvCeent7FYlE5PP5xPg7ewmFQgoG\ngzp58qTq6+vV1dWlnp4eq8vCewqHw3r06JEaGxt18eJFRSIR+mczY2Njb2395/V639oa8JeYugCc\n6F6FSD7RaFSXLl3S0qVLeX7OZjo6OhQMBnX+/HmNjo5qZGRE2dnZ2rFjh9Wl4R3a2tq0cePG+NLg\nli1b1NPTo7Vr11pcGd5Hf3+/ysrKlJGRIWmmf93d3fTPRuaTX0y9g+Xz+d5aswwEAiovLzezBHyA\nyclJXbhwQevXr9f27dutLgcJ2rt3r44dO6b6+nrt2bNHFRUVhCub8Hg8evLkSfy4vb1dhYWFFlaE\nRKxYsULd3d2KRqOSpM7OThUUFFhcFRKxYcMGtba2xnt4586dd+YXU+9grV69Wk+fPlVDQ4OkmcA1\nuzE0kl8gENDAwIDGxsbU3NwsSaqpqVF1dbXFlWE+Zp8lQPKrqqrSs2fP1NDQIJfLpdLSUm3atMnq\nsvCeCgoKVF5erhMnTsjlcmnNmjWqrKy0uiwkIDc3VzU1NTp+/LjS0tJUVFSk3bt3z/k57EUIAABg\nMN7rCwAAYDACFgAAgMEIWAAAAAYjYAEAABiMgAUAAGAwAhYAAIDBCFgAAAAGI2ABAAAY7H/uW7JP\nSVZDxQAAAABJRU5ErkJggg==\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": {}
}
]
}