{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Accuracy of Simpson's Rule" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "import numpy.linalg as la\n", "import matplotlib.pyplot as pt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 0.16666667, 0.66666667, 0.16666667])" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "simpson_weights = np.array([1/6, 4/6, 1/6])\n", "simpson_weights" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is a function and its definite integral:\n", "\n", "$$\\text{int_f}(x)=\\int_0^x f(\\xi)d\\xi$$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "alpha = 4\n", "\n", "def f(x):\n", " return np.cos(alpha*x)\n", "def int_f(x):\n", " return 1/alpha*(np.sin(alpha*x)-np.sin(alpha*0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plotted:" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VGX6//H3MzNpk94oAaQ3KdKLCCIgAkqRtup3LQh2\nxII0WSWogPUnIq7i7tpdG6CggApI1AWDCwJGgUAMHSnpvc08vz8SGGApSSaZM+V+XVcuc5LJnJvb\n5J4zn9OU1hohhBC+w2R0AUIIIVxLBr8QQvgYGfxCCOFjZPALIYSPkcEvhBA+Rga/EEL4GKcHv1Lq\nLaXUcaVU0kUes0gptVcptUMp1dnZdQohhKi+mtjifxsYcqFvKqWGAS201i2Bu4HXa2CdQgghqsnp\nwa+1/hHIvMhDRgDvVjx2MxChlKrr7HqFEEJUjysy/gbAoTOWDwMNXbBeIYQQ5+GqnbvqnGW5ToQQ\nQhjE4oJ1HAEanbHcsOJrZ1FKyYuBEEJUg9b63I3ri3LF4F8JTAY+Vkr1ArK01sfP+8CVl579ubk2\nkpPz2bu3lNRURUpKECUlJjp3LmDIEDvjxllp0SKoZv8FLhYfH098fLzRZbgF6YWD9MJBeuGgVJVm\nPlADg18p9RFwNRCjlDoEzAH8ALTWS7TWq5VSw5RSKUA+MMGZ9YWGmunWLYxu3RxfO3q0iK1bNcuX\nK+bNMxMbm8cNNxQxYUIgXbqEOLM6Q+zfv9/oEtyG9MJBeuEgvXCO04Nfa31zJR4z2dn1XExcXCBx\ncYEMHw5lZZpt24rYuFExYIAf0dG5jBtXzAMPhNKoUUBtliGEEB7BFVGPS1ksiu7dw+jeHWw22LKl\nmO++M/Hqqya6d8/g/vtNjB0bgcmNz1m+4447jC7BbUgvHKQXDtIL5yh3uRGLUkpXJuOvrpwcG19/\nnU1Cgj+geOCBIiZPjiAkxFxr6xRCiNqmlKryzl033u6tWWFhZsaPj+K110K47TYbH35opkGDMh56\nKJ20tFKjyztLQkKC0SW4DemFg/TC4VQvlFI+9VFTvC7quRSloFevMHr1gr178/n0U0XTpprbb09j\nzpxwYmP9jC5RCFEF7pJa1LaaHPw+E/VczL59hfz734Xs3BnCAw/k8MQTkVitEgEJ4e4qYg6jy3CJ\nC/1bqxP1yOA/w969+bz1ViknT/rz1FOF3H13NDX4IiuEqGEy+CXjd1rLlsEsWBDBxIk2nn46kPbt\ns1m7NsfldUiW6yC9cJBeOEgvnCOD/zx69w5l8eJgevWyMW5cAGPHpnHiRInRZQkhPEhycjKdOnUi\nLCyMxYsXG13OWSTquYScnDKWLMnht9+szJ+fzz33SPwjhLtw56hn4sSJRERE8NJLL9XI80nU40Jh\nYRamTYtiypQy5s4NpHfvTJKTi4wuSwjh5g4cOMDll19udBnnJYO/krp2DWHxYiv16tnp1s3Es8+m\nU1sbGpJfOkgvHKQXDu7eiwEDBpCQkMDkyZMJCwsjJSXF6JLOIoO/Cvz9FRMnRjNnThmLFwfQt28G\nhw7J1r8Q4mzfffcdffv25bXXXiMnJ4cWLVoYXdJZZPBXQ+vWVl55JZjgYOjYET78MKNGn79///41\n+nyeTHrhIL1wqGwvlKqZj+py1/0PMviryd9fMXlyFA8+WMqUKUGMG5dGXp7N6LKEEGfQumY+qqsm\nz7atSTL4ndS9eyiLFvmxb5+Jjh0L+PXXfKef093zS1eSXjhILxykF86RwV8DwsIsPPlkFL17l3LV\nVX68/XbNRj9CCM8kUY+XUwpuuimKadNKmDrVyoQJaZSWVu9/umS5DtILB+mFg6f0wl2jHjmBqxZk\nZJTy3HP5WCyKVasC5c5fQtQSdz6Bq6bJCVxuLirKjwULIrjsslK6drWzcWNulX5e8ksH6YWD9MJB\neuEcGfy1xGSCe+6JYezYYoYMCZDcXwjhNiTqcYFff83jxRctTJiQy0svxcq1foSoIRL1SNTjtjp2\nDOH55xVLlwZy/fUZFBbajS5JCOHDZPC7SL16AbzwgpWjRxW9e+de9D6/kl86SC8cpBcO0gvnyOB3\nIavVTHx8JJGRNrp1KyYlpdDokoQQPkgyfgNoDe+9l8H33wexapWNHj1CjC5JCI8kGb9k/B5DKbj9\n9ihuvLGIa6/1Y9WqbKNLEkL4EBn8Bho+PJL77ivi5puDeP/9zNNfl/zSQXrhIL1wkF44x2J0Ab6u\nT59wQkPzmDzZSm5uOvffH210SUIILycZv5vYvbuAefPMzJ6dy7RpMUaXI4RH8OSMPzs7m/79+xMa\nGsratWsJCLj4pV0k4/dCbdpYeeopOwsWhPDEEyeNLkcIUcuee+45Bg4cSFZWFn/88YdL1y2D3400\nbRrEvHmKRYu28sgjaUaX4xYky3WQXjh4ei/sdjvvvPMOd9xxBz/88IPLb8ouGb+badgwgLvvNvHB\nB0GUlKTx2msS+wjhbTZv3ozZbKZ9+/aGrF8yfjd18mQJs2eXMmZMAYsXxxpdjhBuyVMz/vnz5/Pr\nr7/y8ccfV/pnJOP3AbGx/syb58eyZVYefFAyfyG8yffff0/37t0NW78MfjeUlJQAlA//Z57x47PP\nrEyZ4pvD39Oz3JokvXDw5F7YbDZ++uknOnfubFgNMvjdXJ065Vv+n34azMMPyw5fITzd9u3bycvL\n44orrjCsBsn4PcSJE8U8/riN228v4PnnZYevEOCZGf/ixYuZN28ef/75Z5V+TjJ+H1SnTgBPP23i\nrbeCeeYZ2fIXwlNt3rzZsKN5TpHB74ZOZfznql8/kPh4xUsvhbBwYbprizKIJ2e5NU164eDJvUhM\nTDxr8L/88svMnj2bJUuWuKwGGfwepnHjQJ580s6cOSG89Zbcx1cIT5KRkcEff/xBhw4dgPLLNnz6\n6aeMHDmSq666ymV1yAlcbqhDh/4X/X6LFlZmzszj4YdDCA3NYty4CNcUZoD+/fsbXYLbkF44eGov\ntm7dCnB68G/evJlOnTrRo0cPl9Yhg99DtW8fwiOP5HDnnVZCQ3MYMiTM6JKEcDtqbpX2eV6QnlMz\nO5C3bt2KyWSiffv2JCYm8sorrxAXF8fnn3/OjTfeWCPrqAwZ/G4oKSnhklv9AN26hXHffVncdFMQ\n69bl061bcO0X52IJCQkeu3VX06QXDpXtRU0N7Jqybds2WrVqRVBQEL169SIoKIiHH36Ydu3aubQO\nyfg9XN++EYwfX8DQoSa5h68Qbm7Hjh107dr19PKuXbtcfoE2kMHvliqztX+m4cMjufrqQgYOtHH8\neEntFGUQ2cJ1kF44eGIvCgsLSUlJOT34jx8/TkxMDErVTBxVFTL4vcRtt0XRsmUJAwcWkp9vM7oc\nIcQ5du3ahd1uP32Nns2bN9OnTx9DapHB74YudBz/xSgFDz4YRVCQnaFDsykrc69ss7o8+Xjtmia9\ncPDEXvz++++YzWaUUtx7770kJibyl7/8xZBanB78SqkhSqndSqm9SqkZ5/l+f6VUtlJqW8XH35xd\npzg/kwlmzozgxAkzN92UjoedyS6EV9u1axdXXHEFVquVhg0bEhMTY9j1epy6Vo9SygwkA4OAI8B/\ngZu11rvOeEx/4FGt9YhLPJdcq6eG5ObamDGjkFGjinj1Vbmuj/BennStnuHDh9OiRQtefvnlav28\nO12rpweQorXer7UuBT4GRp7nca7fe+HDQkPNxMf78dFHwfy//+cbl3YQwt0lJSUZlumfy9nB3wA4\ndMby4YqvnUkDVyqldiilViulXH/skoepTsZ/rjp1Anj8cTvx8SGsWJHtfFEG8cQst7ZILxw8rRcZ\nGRkcPHjQbQa/sydwVeY91i9AI611gVJqKPAF0Op8D1y48A7q1GkCQHBwBM2adTp9aOOpYSjLVVt+\n8MEu3HZbIC+8sIZWrYJOHwZ36g/H3ZdPcZd6jFzevn27W9Vj5PL27dvxJFu2bKF169bUr1/f6edK\nSEjgnXfeAaBJkybVeg5nM/5eQLzWekjF8izArrV+7iI/sw/oqrXOOOfrkvHXki++yGTVqgC2bDHT\nsGGA0eUIUWM8JeN//PHHycvLY9GiRdV+DnfK+LcALZVSTZRS/sBfgJXnFFVXVZyhoJTqQfmLjVxW\n0oVGjYqkS5ciBg8ukmP8hTDA559/zujRo40u4zSnBr/WugyYDHwD7AQ+0VrvUkrdo5S6p+JhY4Ek\npdR2YCFwkzPr9AU1kfGf6557yo/xHzkyC7u9xp++1nhallubpBcOntCLDRs2cPPNN/Ovf/2LkJAQ\ntzrb2OmLtGmt1wBrzvnakjM+fw14zdn1COeYTDBjRjgzZuRz331pLFkih3kKUZtatWrFli1bSEpK\nYvny5UaXcxa5566PSUsrYfp0G7NmFfDoo9FGlyOEUzwl468J7pTxCw8TE+PPrFl24uNDWbMmx+hy\nhBAGkMHvhmoj4z9Ty5bB3HNPIbfcEkBysntfytkTslxXkV44SC+cI4PfR/XrF8611+YxdGgZ2dll\nRpcjhHAhyfh9mNbw7LMZBATAunWRmM1yZQ3hWSTjl4xfVJFS8OijERw4YGHKlDSjyxFCuIgMfjdU\n2xn/mQICTMyeHchHH4Xy5pvud0E3yXIdpBcOZ/ZCKeUTHzVJbrYuqFPHn+nT85g6NZS2bXPp2zfU\n6JKEqBRfiXlqmmT84rSvv87ks88C2bLFxGWXyTV9hPAEkvELpwwZEkmPHvlcf30RxcUedF0HIUSV\nyOB3Q67M+M91110xaK254w73uI6e5NoO0gsH6YVzZPCLs5hMMGtWCOvXW3n5ZTnSRwhvJBm/OK/k\n5ALmzrXw5ZdF9O8fZnQ5QogLkIxf1JjWra3cfnsB48b5ceRIsdHlCCFqkAx+N2Rkxn+m666LoGvX\nQoYPL6SkxJidvZLlOkgvHKQXzpHBLy7q3nujKCxU3HWX+53cJYSoHsn4xSVlZpby6KNlPPVUAfff\nL9fwF8KdSMYvakVkpB+PPWZnxoxQNm3KNbocIYSTZPC7IXfJ+M/Url0wt9ySz+jRZk6cKHXZeiXL\ndZBeOEgvnCODX1TaDTdE0q5dISNG5GOzSSwnhKeSjF9USWmpZvr0HIYOLeXVV+WG7UIYTTJ+Uev8\n/BQzZwby/vuhfPppltHlCCGqQS7L7IaSkhLo0KG/0WVcUN26AUyenMPddwfRqVMhrVoFVerntNYU\nlRVRUFpAYVkhNruNMnvZ6Q+bLl8GMCszZpOZrZu20rtvb8zKjMVkIcgvCKuflSBLEGaTuTb/mW4n\nISGB/v37G12GW5BeOEcGv6gSrTW5tgzqdDhCh5Ep9H8gh1lPFZBblkVGYcZZH5lFmeSV5FFQWkB+\nST4FpQX4m/0J9g8mwByAn9kPi8mCxWQ5PdhPDXOb3YZN28hNziXgUMDpF4lTLxynnsvqZ8XqZyXY\nP5jIwEiigqKIDIokKrDiv0FRp79eN6QucaFx1A2ui5/Zz+BOCmEcyfjFaVprcsrSOFaSyrHiVNJL\nD5NRepSMkqOklx4p/7z0KAEmK1F+cYRb6nBgdxixwRGMHx5HVFDUWR+RgZGE+IecHs5WP2uNbaVr\nrSm2FZ9+EcgrySOzMJPMoszyF53CzLNegNIL0zmRf4KjuUc5kX+C6KBo4kLjzvpoGNaQZpHNaBbZ\njIZhDbGYZLtIuL/qZPwy+H2MXdtJKznE4eLdHCtO5XhxKn8W/8HximFvVhbqBTSjrn8zYvwbEeUX\nR7Rfg/L/+jcgyq8+ASbr6efLzi7jkUdKiY8v4IEHPOPkLpvddvpF4MyPgzkH2Ze5j9TMVI7nH6dR\nWCOaRTajaUTT0y8IbWLa0Cq6FQEWuVGNcA8y+L1ETWT8dm3nRMl+DhXt5GDh7xws2smhop0cLtqF\n1RxOw8A21PNvTr2AZtQLqPivfzNCLJFVXtfvv+czb54f331XQvfuIU7VfS6jstzismIOZB8gNTP1\n9ItBSmYKu9N2sz9rP43CGnF57OW0jWlb/t/YtrSJaUOIf83++88kubaD9MKhOoNf3st6gSJbPvsL\nfyW1cBupBdtILdzGoaJdhFliaBR4OY0CL6ddSD+GxtxLw8C2hFgianT97doFM25cJqNH+5OUVEZE\nhOf/WgVYAmgV3YpW0a3+53slthJSMlLYdXIXO0/uZE3KGl766SX2pO+hXkg9OtfvTOd65R+d6nUi\nLjSuxm+WLYQzZIvfwxTYcthbsIXUgl9OD/njxftpFNiW5tYuNLN2pllQZy4Lao/V7Nqbpi9YkEFw\nsObbb6PxxTlns9tIyUhh27FtbPtzW/l/j21DoU6/GHSp34WeDXpyWfhl8mIgaoREPV7Gru0cLtrF\n7vxE9uQnkpyfyPGSfTQN6kRza5fyQR/UmYaBbfEz+RtdLoWFdh59tIBJk4qIj5eTu6B8J/SR3COn\nXwh++fMXNh/ZjNaaXg170athL3o26Em3uG6EBrj2hVp4Bxn8Hq7AlsOuvI18/8u/yWx4jL35PxPu\nV4fWwb1OfzQJ6ohFue+hiPv3FzF7tolly4q49lrn79zljVmu1pqD2QfZfGQziYcT2XxkM9uPbadF\nVAt6NuhJn0Z96Ne4H00impz1rsAbe1Fd0gsHyfg9TE5ZGjvz/sPveT/wW+73HClOpqW1B9E0YESd\nh2kd3JMwi2dtOTdpEsiECVnccksAO3aUEBdn/DsRd6OUonFEYxpHNGZ8u/FA+X6DHcd2kHg4kTUp\na5i5fiYWk4V+jfvR77J+9GvcD3fZSBOeT7b4XSi79CS/5n5XPujzviet5BBtQq6kXUg/2oX0o6W1\nG34m7zhM8NVX08nKMpOYGIFJLgxSZVpr/sj8gx8O/HD6I7ckl76X9eXqxlczsNlA2sW2k/0EQqIe\nd1NsL+D3vB/ZkbOO7bnrOF6cSruQfnQIvYZ2If1oZu2EWXnnm66yMs20aTnccEMpCxd61rsWd3U4\n5zA/HviRDfs3sH7fegpKCxjYdCCDmg1iULNBNAxraHSJwgAy+A1m0zb+KNjK9tx17MhZx96Cn2kW\n1JlOYddyReggWgZ3r1Q+7+7X6qmsP/8sYto0xXvvFTFyZHi1nkOyXIdze5Gamcr61PWs27eO9anr\nibHGnH4R6N+kPxGBNXvYrjuR3wsHyfgNkF5yhK05a9ias4ak3A1E+cXRKexaRtWdSruQfi4/pNKd\n1K8fyL33ZjNhQgDbtxdz2WXeEWO5i2aRzWjWtRl3db0Lu7az49gO1qWu4/Utr3Pr57fSLrYdQ1sM\nZVjLYXSN64pJSeYmyskWfxXZtI3k/ES2Zq9mS85qTpYcoHPYdXQNG0qnsGuJ8qtvdIlu54030jl2\nzMzmzeFYLJJJu0JRWREbD25kTcoaVu9dTVpBGkNaDGFYy2EMbj6YqKAoo0sUNUSinlqSU5bG1uyv\n2Zqzmm053xDj34iuYcPoFj6M1sG9vDanryllZZoZM3IZPLiYxYtjjS7HJ+3P2s+avWtYnbKa7/d/\nT4e6HU6/G+hUr5O8G/BgMvhriNaa1MLt/Df7K7bmrOZQ4U46hg6ga/gwuoYNJca/dneieUvGf6YT\nJ0qYOlXz1luFjB5d+exZslyHmupFUVkRPxz4gdV7V7MmZQ05xTkMbzWcEa1HMLDpQIL8Knd/BSPJ\n74WDZPxOKNOl/J77A5uzV5CY9QV+pgB6hA/n/+o/TbuQvl5zmKVR6tTx5777cpg4MYguXYpo0iTQ\n6JJ8VqAlkMHNBzO4+WAWspC96XtZmbySFza9wP8t/z8GNB3AyNYjuaHVDcRY5Ygsb+TTW/yFtjy2\n5XxDYvYXbM1eTf2AFvSIGEmv8FE0Cmwrx0jXgjffzODQITNbtoRJ3u+G0grSWLVnFSv3rGRd6jo6\n1u3IyNYjGdF6xHkvWCeMJ1FPJWSWHue/2V+SmPUFv+f9QJvg3vSMGEXP8BFE+zeo9fX7OpsNZszI\n5pprSnjjDcn73VlRWRHf7fuOFbtXsHLPSiICIxjRagSj2oyiZ8Oesl/ATcjgv4CjRXtJzP6CxKwv\nOFT0O13ChtAzfBRdw4cSbK7e8eW1yRsz/jOV5/12/vGPQsaNu/j1/yXLdTCyF3ZtZ8vRLazYvYIv\nkr8gqyiLG9vcyJi2Y+jbuK/L71YmvxcOkvGf4WDhTjZmfcamzKXklKXRM2IUf6n3BB1Dr5G83mB1\n6vjzwAPlN2vv0qWQ5s3df2eirzMpEz0a9KBHgx7MGziP3Wm7WbZzGY+tfYxD2YcY2XokYy4fw4Cm\nA/A3y/WZ3J3XbPFrrdlfmMSmrKVsylpKoS2XKyPHcmXEWNoE95a3pW7on/9MZ98+C1u3huHnJ3m/\np9qXuY/lu5azbNcydqft5vpW1zOm7Riua36dRxwh5Ol8Luo5ddjlpszyYV+qi7kyYixXRo6llbWH\nDHs3Z7PBrFnZ9O1byptvytEj3uBo7lE+3/U5y3YtY+ufWxncfDBj2o7h+pbXy/0GaolPDH6tNSkF\nW0/HOABXRo6lT8Q4Wli7esWRON6e8Z/p5MnyvP/11wu56ab/zfsly3XwtF6czD/JiuQVLN+1nP8c\n/A/9m/RnTNsxjGg9gsigqt/b+Uye1ovaZEjGr5QaAiwEzMA/tdbPnecxi4ChQAFwh9Z6W1XWobUm\nuWDz6S17i/KnT8Q4ZjRbSrOgTl4x7H1VbKw/kyfncO+9Vrp1K6RFC4kGvEVscCyTukxiUpdJZBVl\n8dWer1i2axkPrnmQPpf1YWzbsYxsM1LOFTCAU1v8SikzkAwMAo4A/wVu1lrvOuMxw4DJWuthSqme\nwCta617nea6ztvjt2s7u/J/KM/vMZQSagukTOY4+keNoHNhehr2XeeutdFJSLPzyi+T93i63OJfV\ne1ezdNdSvv3jW3o06MHYtmO5se2N1AmuY3R5HsflUY9SqjcwR2s9pGJ5JoDW+tkzHvMGsEFr/UnF\n8m7gaq318XOeS3++ooxdeRsrdtAuI9QcTZ+KHbSXBV1e7TqF+7Pby/P+K68s5Z//lC1AX5Ffks/X\nKV+zdNdS1uxdQ+f6nRnbdiyj246mfqhc8LAyjBj8Y4HrtNZ3VSz/FeiptX7wjMd8CSzQWm+qWF4H\nzNBabz3nuXT407GEm+vSPfgGultvoL5/i2rX5sl27dpI27Z9jC7D5TIzy4iPj2LevGOMGWMFYOPG\njfTp43u9OB9v70VhWSHfH/qer1K/Yu2BtbSJasMNzW5gWLNhNAg5++RKb+9FZfn5+REdHe3yjL+y\nrxrnFnXen+uw4XKaNWwJpHEgeCl+LVrRqVNXALZvL3+d8IXltLTjpKUtd5t6XLk8dWoXZs26gv/+\n9y0iI8t/bU6ePEBycvn3W7cuf7wvLh86tIeTJxu6TT21tTyItsSd7MfBI7tYnZ3Isz+9QOjBKFr5\nd2Fwu1uJscSxbt2P/Oc/x92iXlcvJydv5aefvgIgNrZ6J6A6u8XfC4g/I+qZBdjP3MFbEfUkaK0/\nrli+YNTjLkcYCWM99lg6n3/uxwsvhEreLyjTpSTlbmBj5lISsz+njn8T+kSMpXfEGOICfTMVOOXw\n4a+4//7hLo96LJTv3B0IHAV+5uI7d3sBCy+0c1cGv4DyvP/KK7MIDy9j8mTJ+4WDTZfxW94PbMpc\nyk9Zy4n0q0+finN3Gga2Nro8l6vu4HfqDCetdRkwGfgG2Al8orXepZS6Ryl1T8VjVgOpSqkUYAlw\nvzPr9AUJCQlGl2AokwmWL7eydWsI7733pdHluI2kpASjSzCcWVm4InQAV2WP5+0OR7ir4StklR3n\nb3sH8ODODnz051wOFv6ObERenNPH8Wut1wBrzvnaknOWJzu7HuFb4uL8ef/9HMaMCWDAgCIaNpTr\n94uzmZWZ9qH9aB/aj0kNF5Kcn8imrKXEpwwl0BRccWLnWJoEdZTDv8/hVmfuukstwn3MnJnGJ58E\n8OKLIfj7yx+vuDStNXsL/svGzKVszPoMs7JwZcRY+kSOpXlQF696ETAk469JMvjF+WgNV12VSWCg\nnYcfjja6HOFhyq/nte30i4BNl53eJ9DK2sPjXwQMyfhF7fD1jP9M33+fwPLlwSQlWVm7Nsvocgwl\nGb9DZXuhlKK5tQu3NZjPG5fvYXazFfibgnhl/x1M/K0x/zz8CDvzNmLX9tot2M147fX4hfeoW9ef\njz7K4cYbrbRqVUTjxpL3i6pTStHU2pGm1o78X9xTHCzcyaaspbx+8D5ybWn0jhhDn4ixtA25CrMy\nG11urZKoR3iMJ544ybvvBvHSS1YCA+XNqqg5h4uS2ZS1jE2ZS8koPUqviFH0DB9Jh9Br8De574aG\nZPzC62kN11yTgdbw2GNRRpcjvNSfxX/wU9Zyfs5ayf7CX+kYOpCeESPoFnY94X7udZ9oyfi9iGT8\nDmf2QilYujSUPXsCWb3a9/J+yfgdarMX9QOaM7ruNJ5t/SNvtv+DXhGj+G/2V9y7syXTk/uw9Niz\nHCzc6dHnCkjGLzxKTIwfn3xSyPXXW2nTppBmzeT6/aL2hFliGBB9GwOib6PUXkxSXgI/Z60kPmUI\nFpM/PcKH0yN8BJeHXIVF+RldbqVJ1CM80rx5abz2WhALFwYRFCRvXIVraa3ZV7iDn7O/5OfslfxZ\nnMIVoQPpEjaULmHXEePf0CV1SMYvfIrWMGRIOrm5JmbNcu42fkI4K6P0T7blfMsvOV+zPedbovzi\nyl8EwodweXAf/EwBtbJeyfi9iGT8DhfqhVLw8cdhHDzoz4oVvpH3S8bv4G69iPKrz8Do25nW9CPe\n63iCBxr/gwCTlfePzOLWX+vwzB8jWX3ydY4V7zO6VEAyfuHBIiP9WLq0iGuvtdK2bQGtWlmNLkkI\nzMpMm+BetAnuxS1x8eSUpbEtZy2/5HzNx3/OJcBkpWPoADqGDqRj6DVE+tVzeY0S9QiP9+KLJ3nx\nxWBefjmAkBDvPvFGeDatNYeKdvJr7nfsyF3Pb3nfE+UXxxWhA+kYOoD2If0JsURU+vkk4xc+bcSI\nNI4dM/O3v0Xi4ZdfET7Epm2kFvxy+oUgOf8nGgS2oWPoANqF9KVN8JWEWi58zopk/F5EMn6Hyvbi\nww8jOHGbeDhsAAATv0lEQVTCj2XLMmu3IAO5W65tJG/phVmZaRncnTH1ZvBUy2/5oGMadzZ4EX8V\nxMoTC5n0W2Me3NmBvx+8j4SMDzlRfKBG1isZv/AKoaEWli+Hq68Opm3bfNq1Cza6JCGqzM8UQPvQ\nq2kfejVQftvJfQU72Jn/HxKzPuetw1OxKH8uD7mKNsFXEllaWq31SNQjvMprr6UTHx/EwoX+hIXJ\ndo3wLlpr/iz+g535P7I77yd2Z2/g4N9SJOMXYvz4NPbuNTN3ruT9wrtJxu9FJON3qE4v3n03irw8\nCx99lF7zBRnIW3LtmiC9cI4MfuF1goJMfPGFhVWrQvnllxyjyxHC7UjUI7zW++9nMmVKEC+9ZCI2\n1t/ocoSocRL1CHGOW2+N5C9/yWXBgiLKymSjQohTZPC7Icn4HZztxeLFMYSH23njDc8/vl9ybQfp\nhXNk8AuvZrEoVq4M4pdfgli7NtvocoRwC5LxC5+wfn02o0YF8cwzNrl5i/AakvELcREDB4bz+OPZ\nPPusnbw8m9HlCGEoGfxuSDJ+h5rsxaxZsXTvXsiLL+bgiW8uJdd2kF44Rwa/8CkffRRJdraFDz/M\nMLoUIQwjGb/wOTt3FtC7t5mHHiqhe/dQo8sRotok4xeiki6/3Mobb+SzaJE/x48XG12OEC4ng98N\nScbvUFu9uPnmKG69NZf584spKfGMd5qSaztIL5wjg1/4rIULo6lb18Zrr0neL3yLZPzCp504UcoV\nV5QyfHgRw4df+BZ3QrgjyfiFqIY6dfxYsULz0Uch7NiRZ3Q5QriEDH43JBm/gyt60aNHMK+8ksdL\nL1k4caKk1tdXXZJrO0gvnCODXwjgzjujuOWWXObN85ydvUJUl2T8QlSw2TTXXJOJ3a6YPl1u2yjc\nn2T8QjjJbFZ88UUoBw/68emnnn8ZZyEuRAa/G5KM38HVvYiK8uOrrxRffmnl559zXbruS5Fc20F6\n4RwZ/EKco2PHYP75zwIWLfLn8GE5s1d4H8n4hbiAWbNO8v77Qbz4YhDBwWajyxHif0jGL0QNmz8/\nlo4di3nhBc+8jLMQFyKD3w1Jxu9gZC+Ugs8+iyAvz8Lbb6cbVscpkms7SC+cI4NfiIsIDjazZo0f\nP/4YzDffyJE+wjtIxi9EJfzwQy7XXx/ArFkldOgQYnQ5QgCS8QtRq/r1C2XRojyef97C0aNypI/w\nbDL43ZBk/A7u1IsJE6K4++5cnnqq1JAbtkuu7SC9cE61B79SKkoptVYptUcp9a1SKuICj9uvlPpV\nKbVNKfVz9UsVwnjPPRdL165FLFiQi831s1+IGlHtjF8p9TyQprV+Xik1A4jUWs88z+P2AV211he9\n24Vk/MJTFBXZ6dUrl9hYG1OmyDX8hXGMyPhHAO9WfP4uMOoij5XLXQmvERho4uuvg9i5M4ClS40/\nzFOIqnJm8NfVWh+v+Pw4UPcCj9PAOqXUFqXUXU6sz2e4U65tNHftRb16/qxerVixIoT//CfLJeuU\nXNtBeuEcy8W+qZRaC9Q7z7dmn7mgtdZKqQvlNH201n8qpWKBtUqp3VrrH8/3wDvuuIMmTZoAEBER\nQadOnejfvz/gGACy7FvLp7hLPecu//vfXbjpJiuZmWto0iSIDh3Kv39qMNXkcmrq9lp9fk9aTk3d\n7lb1uHI5KSmB9evfASAwsHo3DnIm498N9NdaH1NK1Qc2aK3bXOJn5gB5WuuXzvM9yfiFR3rzzXRm\nzAhmwQJo0CDQ6HKEDzEi418J3F7x+e3AF+c+QCllVUqFVnweDAwGkpxYpxBu5+67o7n//lzi421k\nZZUZXY4Ql+TM4H8WuFYptQcYULGMUipOKbWq4jH1gB+VUtuBzcBXWutvnSnYF7hrrm0ET+nFvHmx\nXHttAXPnFlBUZK+VdUiu7SC9cM5FM/6LqTg8c9B5vn4UuL7i81SgU7WrE8KDvPNODIMHZ7BgQTZz\n5kRiktMjhZuSa/UIUYMKCmz06pVHbGwZDz8cbXQ5wsvJtXqEcANWq5n1662kpATywQdyjL9wTzL4\n3ZCn5Nqu4Im9iI31Y906M+vXB7NyZc1dyllybQfphXNk8AtRC1q2DOSbb2x89pmVdetcc4KXEJUl\nGb8QtSghIYcRIwKZPLmI3r3DjC5HeBnJ+IVwQ/37h/HBBwUsXhzIjh35RpcjBCCD3y15Yq5dW7yh\nFyNGRPD3v5ffxGXPnoJqP4/k2g7SC+fI4BfCBf761yjmz8/l6adN7N9fZHQ5wsdJxi+EC82bl8bL\nLwczf76ifn25ro9wjmT8QniA2bNjmDAhnyeftHPihNy7VxhDBr8b8oZcu6Z4Yy9eeCGG8eMLmD27\njJMnK39ZXcm1HaQXzpHBL4QBFi2KYcyYAmbPLiUtrXrXVBe+TWvYvv1C97+6OMn4hTDQvfeeZOVK\nK/Pn+xMd7Wd0OcJDaA2vvprOvn3+/PFHmGT8QniS11+P5YYbCpk9u5j0dNnyF5emNbzySjr79/uz\naVP1DhCQwe+GvDHXri5v74VSsGRJDEOHFjF7dslFYx/JtR18tRdaw8KF6Rw+XD7069Sp3rtEGfxC\nGEwp+Mc/Yhg1qoiZM8s4dkyO9hH/y2aDF19M5+hRPzZuDCQ2tvrRoGT8QriRadNO8vbbIcydq7js\nMjnOX5QrLdU8+2wWxcUmEhKCiYx03ENLKSUZvxCe7IUXYpk8OZ8nnoDU1Opf3kF4j6IiO3PnZmMy\nwcaNIWcN/eqSwe+GvD3Xrgpf7EV8fAyzZuUzZ46ZXbscF3bz1Vz7fHylF/n5Np58MpfoaDsJCeGE\nhJhr5Hll8AvhhqZNi2bBgnyeecbCli25RpcjDJCdXcrs2fm0aFHGt99GEhhYc+NaMn4h3Ninn2Yx\ncaKVO+8sYNCgCKPLES5y9GgRc+eWMWBAMe+9F43pIjO/Ohm/82GREKLWjB8fQWxsDqNHB5KZmcG4\ncVFGlyRqWXJyAfPnm5g0qYjnnotBVWmkV45EPW7IF3PtC5FewDXXhPH99za++upnlixJR94Ye2/G\nv3lzLnPnWoiPz+f552tn6IMMfiE8QseOwbzxhonkZH+eey6dkhKZ/t5m9epMFi0K4N13C5gyJbpW\n1yUZvxAeJCurjGHDcklLM/HEE8FEREha6+m0hnffTSchIZiVK0vp0ye0Sj8vx/EL4eUiIiz88EME\nvXqVMHVqCamphUaXJJxQVGTnmWcy2bHDn8REXeWhX10y+N2Q5NoO0guHU72wWBTvvRfL1KkFPPGE\niY0bc4wtzADekPGfOFHCtGn5hIRotm2z0rJlkMvWLYNfCA81Y0YM779fyBtvBPDxx7LT15MkJeXx\n2GN2hg8vYv36KEJDa+bErMqSjF8ID5eUVMANN9ipX7+ERx+NIChItufc2VdfZfDvf4ewcGEeEyc6\nf3hudTJ+GfxCeIGcHBvjx2eRlBTIzJlmmjSRC7y5m6IiO6++msXevYEsW2bnyitDauR5Zeeul5Bc\n20F64XCxXoSFmVmzJpq77y5g9mwT69dnua4wA3haxn/wYBFTpxZgMsFvv/nV2NCvLhn8QngJpWDO\nnFiWLSvkww8DePVVOd7fHaxbl8XMmSZuv72QDRui3OIWmxL1COGFjh4tYezYfA4ftjB1qoVmzVx3\nxIgoV1hoZ8mSTH791cqHH5YweHB4raxHoh4hBABxcf5s3BjJpElF/O1vZj77LFOO+nGh33/P56GH\nCgFFUpK51oZ+dcngd0OSaztILxyq2gul4MknY9mwoYSNGy3Mnp110Xv6ehJ3zfhtNnjrrQzmz/dj\n1qxC1q2Lol49f6PL+h8y+IXwct26hfD778H06FHGww9r1qzJkq3/WrBvXwGPPprDgQMmtmyxMWVK\n7V1kzVmS8QvhQ775Joe77jIRHl7G5MlBxMUFGF2Sxysp0XzwQQbr1oXy0EM5zJ0bjdnsuokvx/EL\nIS6psNDG9OkZvP12OCNG5DF+fBRm15446jW2bs1jyRJF48YlvPNOIG3auH4nuuzc9RKSaztILxxq\nqhdBQWZefTWWhIQSkpLMPPxwLjt25NXIc7uK0Rl/ZmYpzz6bwaJFFubOLeKnnyINGfrVJdd0FcJH\ndesWwo4d8PLL6cybF0yLFplMmmSlfn2Jfy6kpETz6aeZrFkTwrBhdlauNBMdXbvXzq8NEvUIIcjM\nLGXmzCw+/DCcq6/O5a9/jSAkRPKfU7SGr7/O4tNPA2jZspCFC/3p2dPYs29PkYxfCOGUPXsKmTKl\ngMTEYIYNy+fGGyMJDPTtRDgxMZsPPzRhscDzz5cxenSk0SWdRTJ+LyG5toP0wsEVvWjVKoivv45m\n5cpSDhwwcffdJXz8cTpFRfZaX3dV1HbGrzVs2pTDI4/k8K9/+TF1agnJySFuN/SrSzJ+IcT/6Ncv\nlI0bYcOGXB5/3MykSaVcc00eo0eHERlp/LVmaktZmea777JYtcpMSYmZadOKuP/+UPz8rEaXVqMk\n6hFCXFJiYh5PP11MQkIY3brlMHx4EK1be88wzM4uZeXKHL77zkqdOqU89FAZEydGYrG46RlYZ5CM\nXwhRq1JTi3jhhTw++SSEyMgSBgwoY+DAcIKDPW9HsNbw8885rF1bRlJSKFddlcv06RYGDQozurQq\ncWnGr5Qap5T6XSllU0p1ucjjhiildiul9iqlZlR3fb5Ecm0H6YWDO/SiWbNAXn89hmPH/Pnb32wk\nJSnuvNPO009nkJCQ5bLLQFc349e6/AJqf/97OhMnFvLee2YGDbKTkqL59tsojxv61eXMzt0k4Ebg\nhws9QCllBhYDQ4DLgZuVUm2dWKdP2L59u9EluA3phYM79cLf38SECZEkJkaSnGzn+uvtrF5t4rbb\nypg7N5Mvv8wgI6P2LgiXmlr5XhQX29m4MYdXXkln0qQCXn7ZRP36mqVLyzhwIJh582KIi3O/C6nV\npmrv3NVa74bytxkX0QNI0Vrvr3jsx8BIYFd11+sLsrK8++5JVSG9cHDXXjRqFMCcOQHMmQP79hXx\n2Wd2Vq408f77JmJi8mjZspj27U106GClbt2AGrlwWX7+hXuRn28jKSmf334rJSXFTGpqCE2awLXX\nap54wka/fqEo5Tln2daG2j6qpwFw6Izlw0DPWl6nEMIgTZsGMn16INOnl29pJyTYWLtW8+OPinfe\nUUAxl11WSP36NuLiFA0amGnUKIA6dQLw96/aK8LJkyUcPVrC0aOlHDli4/BhM0eOBJCREUDz5tCj\nh+axxzRDhtipX983IpzKuujgV0qtBeqd51uPa62/rMTzy97aati/f7/RJbgN6YWDp/UiIMDEddeF\nc9115ctaQ0pKIRs3apKSNHv2wC+/wJEjkJOj8fcvIzS0lMBAG/7+GotFYzZrbDZV8QEFBWby8vzI\ny9vHqlVQp46mYUNN8+YwaJCdbt1sdOmiCAiQQX8xTh/Vo5TaAEzVWv9ynu/1AuK11kMqlmcBdq31\nc+d5rLxICCFENVT1qJ6ainoutNItQEulVBPgKPAX4ObzPbCqhQshhKgeZw7nvFEpdQjoBaxSSq2p\n+HqcUmoVgNa6DJgMfAPsBD7RWsuOXSGEMJDbnMAlhBDCNVx6kbbKnMyllFpU8f0dSqnOrqzPlS7V\nC6XU/1X04Fel1EalVEcj6nSFyp7kp5TqrpQqU0qNdmV9rlTJv5H+SqltSqnflFIJLi7RZSrxNxKj\nlPpaKbW9ohd3GFBmrVNKvaWUOq6USrrIY6o2N7XWLvkAzEAK0ATwA7YDbc95zDBgdcXnPYFEV9Xn\nyo9K9qI3EF7x+RBf7sUZj/sO+AoYY3TdBv5eRAC/Aw0rlmOMrtvAXsQDC071AUgHLEbXXgu96At0\nBpIu8P0qz01XbvGfPplLa10KnDqZ60wjgHcBtNabgQilVF0X1ugql+yF1vonrXV2xeJmoKGLa3SV\nyvxeADwILAVOurI4F6tML24BlmmtDwNordNcXKOrVKYXfwKnjtsMA9J1+X5Fr6K1/hHIvMhDqjw3\nXTn4z3cyV4NKPMYbB15lenGmicDqWq3IOJfshVKqAeV/9K9XfMlbd0xV5veiJRCllNqglNqilLrV\nZdW5VmV68Q+gnVLqKLADeMhFtbmbKs9NV16Pv7J/rOce1umNf+SV/jcppa4B7gT61F45hqpMLxYC\nM7XWWpVfI8RbD/2tTC/8gC7AQMAK/KSUStRa763VylyvMr14HNiute6vlGoOrFVKXaG1zq3l2txR\nleamKwf/EaDRGcuNKH9luthjGlZ8zdtUphdU7ND9BzBEa32xt3qerDK96Ap8XHFdqBhgqFKqVGu9\n0jUlukxlenEISNNaFwKFSqkfgCsAbxv8lenFlcA8AK31H0qpfUBrys8f8iVVnpuujHpOn8yllPKn\n/GSuc/9wVwK3wemzfrO01sddWKOrXLIXSqnLgOXAX7XWKQbU6CqX7IXWupnWuqnWuinlOf99Xjj0\noXJ/IyuAq5RSZqWUlfKdeTtdXKcrVKYXu4FBABWZdmsg1aVVuocqz02XbfFrrcuUUqdO5jID/9Ja\n71JK3VPx/SVa69VKqWFKqRQgH5jgqvpcqTK9AJ4EIoHXK7Z0S7XWPYyqubZUshc+oZJ/I7uVUl8D\nvwJ24B9aa68b/JX8vZgPvK2U2kH5Rux0rXWGYUXXEqXUR8DVQEzFSbNzKI/8qj035QQuIYTwMS49\ngUsIIYTxZPALIYSPkcEvhBA+Rga/EEL4GBn8QgjhY2TwCyGEj5HBL4QQPkYGvxBC+Jj/Dzxlqa49\nxAjcAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_x = np.linspace(0, 1, 200)\n", "\n", "pt.plot(plot_x, f(plot_x), label=\"f\")\n", "pt.fill_between(plot_x, 0*plot_x, f(plot_x),alpha=0.3)\n", "pt.plot(plot_x, int_f(plot_x), label=\"$\\int f$\")\n", "pt.grid()\n", "pt.legend(loc=\"best\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This here plots the function, the interpolant, and the area under the interpolant:" ] }, { "cell_type": "code", "execution_count": 61, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8TfcbwPHPNwkhkRCCEBSxqxpas7RRs6gqVaVm+2tt\nWlvVKDVKtarV6qZ20apVW2pTe+8VewZRROT7++OEayQkueOc5D7v1+u+5Oaee86TR/Lk5Dnf8/0q\nrTVCCCHch4fZAQghhHAtKfxCCOFmpPALIYSbkcIvhBBuRgq/EEK4GSn8QgjhZuwu/EqpX5RSZ5VS\nOx6zzWil1AGl1DalVEl7jymEECL5HHHG/ytQM6EXlVK1gAJa64LA+8B3DjimEEKIZLK78GutVwKX\nH7NJXWB83LbrgUxKqez2HlcIIUTyuKLHHwxE3Pf8BJDLBccVQggRD1dd3FUPPZd5IoQQwiReLjjG\nSSD3fc9zxX3uAUop+WUghBDJoLV++OT6sVxxxj8baA6glCoHRGqtz8a75QDjUaNVDbTW8T6uXNGs\nXKn5+mtNs2aavHk1WbNq6tfXfP+95sSJ+N+Xkh79+/c3PQarPCQXkgvJxeMfyWH3Gb9SagrwEhCo\nlIoA+gNpALTW32ut5yulaimlDgLXgVaP25/fSj869uyY4Ov+/lCxovG4KyICwsPh77+hd2/InRvq\n14dGjaBwYXu/Qtc7evSo2SFYhuTCRnJhI7mwj92FX2vdOBHbdEjMvorvLM7RkKNcyHYhSTHkzg3N\nmhmPmBhYuxZmzoTKlSEoCJo0gebNIVu2JO1WCCFSJUvduTvk0yF0a9yNjvM7svvc7mTtw8sLKlWC\nUaOMvwS++AJ27zbO/Bs2hEWLIDbWwYE7WMuWLc0OwTIkFzaSCxvJhX1UcntEjqaU0rP3zgZg0o5J\nbD69me1ttpMuTTqH7P/KFZg8Gb7/HqKj4cMPoWlTSJ/eIbsXQghTKKXQFry4m2RNijfBN40vLWa1\ncNg+M2aEtm1hyxYYMwb++gvy5oVBg4xfClYSHh5udgiWIbmwcWQulFLySIEPR7Fk4VdK0bV8V5Yd\nXcbYjWMdvG+j9z93rnFB+MABKFAABg6EyEiHHkoISzN7NIo8nD96JyGWbPXctfPcToauGso/Lf+h\nZA7nze22fz8MHgzz5xujgtq3B29vpx1OCNPFtQfMDkMkQUL/Z6mm1XNX8WzFeb3I67w+7XWu3brm\ntOMUKgTjx8M//8DSpVCsGMyYAfJzIYRIjSxd+AEaFG1ANt9sNJ7Z2OlnKMWKwbx5xgXgQYOM0UHr\n1zv1kPGSvraN5MJGciEcxfKFXylFl3Jd2Hx6M0NXDXXJMatWhc2b4d134fXXjYvC0v8XwnX27dtH\naGgo/v7+fPPNN2aHk+pYusd/v0OXDtEvvB9/NvqTKvmruCyuyEjo1QvmzDHuDXjjDeMCsRApmdV7\n/O+++y6ZMmVi5MiRZodiGW7T479fSOYQ3in5Do1mNOLE1RMuO26mTDB2LEybBv37w6uvwrFjLju8\nEG7p2LFjFCtWzOwwUq0UU/gBquSrQung0tSdUpfbd2679NgVKxr3AJQtC889Bz/95LyLv9LLtZFc\n2LhLLl5++WXCw8Pp0KED/v7+HDx40OyQUp0UVfgB3i/5Pv/d/o+289q6/Nje3tC3rzH+/+uvjf7/\n+fMuD0OIVG3ZsmVUqlSJMWPGcPXqVQoUKGB2SKlOiiv8Xp5e9HqhF3/u/ZNftvxiSgzFi8OGDcYw\n0GefNcb/O1JYWJhjd5iCSS5sXJ0LpRzzSC4rX4NI6VJc4QfI4pOFLuW68MGCD9h8erMpMXh7w/Dh\nMGUKtGtnPG7cMCUUIZxCa8c8ksuRUxSIB6XIwg9QInsJ6hepT72p9bhy07zJdl56CbZtg0uXoEIF\nOHTI/n26Sy83MSQXNpIL4SgptvADNCjWgNz+uak/rb6pfxZmzGic+b/7LpQvb0wAJ4Swj7R6nCdF\nF36AzmU7c+DSAfos62NqHEpBhw4wezZ06gQ9exqLwiSH9LVtJBc27pYLafU4T4q5getxTlw5Qa9l\nvRj32jheL/q6gyNLugsXjLn+b9yA6dNl5S9hPVa/gUs8yi1v4HqcXBlz0b50e1r91Yo95/eYHQ6B\ngcZInxdfNMb9b9+etPdLL9dGcmEjuRCOkioKP0D5XOV5peAr1J5cm6s3r5odDh4exkRvQ4dClSrS\n9xdCWEeqaPXcpbVmyKohZEibgcXNFlumR7hhA9Svb1wD6NlT5voR5pNWT8ojrZ4E3J3J8+Clg/RY\n3MPscO4pU8aY3nnGDGjRAm7dMjsiIYQ7S1WFHyB9mvT0fqE3P27+kSk7p5gdzj3BwbBiBVy/DjVr\nPn6dX+nl2kgubCQXwlFSXeEHyOmfk85lO9N6Tmu2ndlmdjj3+PjA778bUz5UqgQnT5odkRDCHaWq\nHv/Dpu+ezvKjy9neZjsB6QMcum97aA0jRsCYMfD338bKX0K4kvT4Ux7p8SdSw2INyZcpH69OeZVY\nHWt2OPcoBT16wKefQuXKsGqV2REJIdxJqi78AJ3KdOJM1Bk6zO9gdiiPaNYMJk40RvzMmWP7vPRy\nbSQXNpILawkLC+Pnn382O4xkSfWF39vLm48qfsTUnVP5adNPZofziGrVjAXe33vPWOVLCHeXN29e\nli5d+sTtzC68SimHDBkPDw8nd+7cDogo8bxcejSTZPXNSvfy3flw0YcUCizEi0+9aHZIDyhdGhYt\nMkb7/PcftGoVZnZIluFu89M8jrvkIrEF1d6iGxsbi4dHqj/3jZfbfNXFsxenxbMtqD+tPscirbdo\nbokSxspe/fsbq3sJ4c601owbN46KFSvSvXt3MmfOTP78+VmwYAEAffr0YeXKlXTo0AE/Pz86deoE\nwN69e6lWrRpZsmShSJEiTJ8+/d4+W7ZsSdu2balVqxYZMmRg+fLltGzZkjZt2lC9enX8/f0JCwvj\n+PHj996zZs0aSpcuTaZMmShTpgxr166NN95Dhw7x8ssvExgYSNasWWnatClX7huznTdvXkaOHMmz\nzz5LpkyZeOutt7h16xbXr1/nlVde4dSpU/j5+eHv78+ZM2eckdIHuE3hB6gRUoMXcr9AjYk1uB59\n3exwHlGokDHWf9iwcIYMMTsaa5C+to075eLu2fyGDRsoUqQIFy9epEePHrz77rsADB48+N7yjNeu\nXWP06NFcv36datWq0bRpU86fP8/UqVNp164de/bY5u+aMmUKffv2JSoqiooVKwIwefJk+vXrx4UL\nFwgNDeXtt98G4NKlS9SuXZsPPviAS5cu0aVLF2rXrs3ly5fjjblPnz6cPn2aPXv2EBERwYABAx74\neqZPn87ChQs5cuQI27dvZ9y4cfj6+rJgwQJy5szJtWvXuHr1KkFBQc5I6QPcotVzv3dLvsvAFQOp\nP60+C5ousMy0DnflzQujR0O/fhAdDfd97wjhUuoTx/xs6P7JHzb61FNP3Sv2zZs3p127dpw7d45s\ncVPe3j+8ce7cueTLl48WLVoAEBoaSv369Zk+fTr9+vUDoF69epQvXx4Ab29vAOrUqXPvl8DgwYPJ\nmDEjJ06cYPny5RQuXPjeL4K33nqL0aNHM3v27HvHuCskJISQkBAAAgMD+fDDDxk4cOAD23Tq1Ole\nUX/11VfZunXrI1+Dq7hd4VdK0b18d3ou7cmHCz9kVM1RZof0iAYNwqhYEV5+2ZjsLe571i25S187\nMVydC3sKtqPcf/br4+MDQFRU1L3Cf/+J27Fjx1i/fj0BAbZ7dmJiYmjevPm9bXPlyvXA/h/+nK+v\nL5kzZ+bUqVOcPn2aPHnyPLD9U089xalTpx6J8+zZs3Tu3JlVq1Zx7do1YmNjyZw5c4JfS/r06ePd\nj6u4VavnLp+0PvSp2Ifftv3GT5utN9IHIHt2WLYMpk41ZvkUQjzo4b/W8+TJw0svvcTly5fvPa5d\nu8aYMWMS3IfWmoiIiHvPo6KiuHTpEsHBweTMmZNjxx68Hnjs2DGCg4Mf2c9HH32Ep6cnO3fu5MqV\nK0yYMIHY2MTdO2RG18EtCz9AkF8Q3St058OFH/LPsX/MDucBd3u5d4v/5MnGzV7uyJ362k/iTrlI\nTPsje/bsHLpvkes6deqwf/9+Jk6cyO3bt7l9+zb//vsve/fufew+58+fz+rVq4mOjqZv376UL1+e\n4OBgXnnlFfbv38+UKVOIiYlh2rRp7N27lzp16jyyj6ioKHx9ffH39+fkyZOMGDEi0V9r9uzZuXjx\nIlevum46ebct/ADFsxWn5bMtaTCtAUcuHzE7nHgFBRnFf+JEGDzY7GiEcI27QzofPhu+/3nnzp2Z\nMWMGmTNn5oMPPiBDhgwsWrSIqVOnEhwcTI4cOejduzfR0dEP7PPh/TVp0oRPPvmELFmysGXLFiZO\nnAhAlixZmDt3LiNHjiQwMJDPP/+cuXPnPtLCAejfvz+bN28mY8aMvPrqqzRo0OCxZ/L3x1KkSBEa\nN25M/vz5yZw5s0tG9aTquXoS6+ctP7Pnwh42vbeJDN4ZTInhSU6fhrAw+N//oHt3s6MRKZ3M1WNo\n1aoVuXLlYlAK6KfKXD0O9k7oO2T0zki9afUsNafP/XLkgKVL4dtvYexYs6MRInVw119+UvgxfmN2\nLdeVw5cP03ZuW7PDSbCXmysXLFli9Pvj/hpN9dypr/0kkgvHc9S0CymN2w3nTIhPWh/6vtiXHot7\nUChLIbpW6Gp2SPEKCTGmd6hSBTJkgHr1zI5IiJTr119/NTsEU0iP/yGHLh2iX3g/xtUbR/2i9c0O\nJ0GbNsErr8CkScZEb0IkhfT4Ux7p8TtRSOYQOpbpSMtZLfn35L9mh5Og556DP/6AJk1gzRqzoxFC\npCRS+ONRLlc5Gj3diNqTa3P8yvEnv8HBEtvLrVgRJkww5vO/bzqSVEX62jaSC+Eo0uNPQN3CdTkT\ndYZqv1VjU+tNZEhrzWGeNWsayzjWrGmc+cdzU6EQ8XLHi5rCID3+x9BaM2TVELw9vVnWYhmeHp5m\nh5Sg4cONs/8VKyDAOssLCyGcTHr8DqaUolv5bpyOOk2rv1qZHc5jde8OVavCa6/BzZtmRyOEsDIp\n/E/g7eXNx5U+ZtGhRQwMH/jkNzhAcnq5SsHIkUar5+234c4dx8dlBulr20gubCQX9rG78Culaiql\n9iqlDiilesbzephS6opSakvc42N7j+lqAekD6FupLyPXjWTCtglmh5MgDw8YNw6uXIEOHcAiXTwh\nhMXY1eNXSnkC+4CqwEngX6Cx1nrPfduEAV201nWfsC/L9fgftvXMVoavGc4fb/5BtRDrDp6/etWY\n16dePfeey18Id2BGj78McFBrfVRrfRuYCrwWX2x2HscSQoNCea/UezSc3pCtZ7aaHU6C/P1h/nzj\n7H/8eLOjEUJYjb2FPxiIuO/5ibjP3U8DFZRS25RS85VSxew8pqkq561M/aL1qTmxptMWbXdE/zIo\nCObNgx49jEXcUyrp5dpILmwkF/axdxx/YvpEm4HcWuv/lFKvALOAQvFtOKrXKLIFG0uq+fr5kr9o\nfp4p+wwAO9bvALDE89eLvM6+jfso37c8u4bvIiB9wL1vxLvL41nl+ZQpYTRqBCNGhJMnj/nxJPX5\nXVaJx8znW7dutVQ8Zj6/u16tVeJx5fPw8HDGjRsHQN68eUkOe3v85YABWuuacc97A7Fa688e854j\nwHNa60sPfd7yPf77aa0ZsWYEN2JusPqd1Xh7eZsdUoJ++QWGDIG1ayFrVrOjEUI4khk9/o1AQaVU\nXqVUWqAR8ED1VkplV3G3CCqlymD8srn06K5SFqUUXcp1IfpOtKXn8Qd45x14803jYq+M8RdC2FX4\ntdYxQAdgIbAbmKa13qOUaq2Uah232RvADqXUVmAU8JY9x7QSL08v+lTsw74L+3hv9nsO2+/DbQ5H\n+PRTYz7/Vq0gkWtAW4IzcpFSSS5sJBf2sXscv9b6b611Ya11Aa310LjPfa+1/j7u4zFa6+Ja61Ct\ndQWt9Tp7j2klPml96P9if+YdmEf/5f3NDidBHnFj/I8ehf7WDVMI4QIyV4+DRFyJoM/yPnxW9TPe\nf+59s8NJ0LlzUL68Mb6/RQuzoxFC2Evm6jFR7oy56fVCL7ot6sbM3TPNDidB2bLB3LnG3D6rV5sd\njRDCDFL4HahY1mJ0LtuZVn+1Yunhpcnej7P7l0WLGjd2NWwIx12/3ECSSC/XRnJhI7mwjxR+BysT\nXIb/lfwfDX5vYOkVvF55Bbp2NWbzvH7d7GiEEK4kPX4nmbN/Dn/s+YPV76ymcGBhs8OJl9bGUM+o\nKJg2zbgALIRIWaTHbyGvFnqVaiHVqDy+MhFXIp78BhMoBWPHwsmTMGiQ2dEIIVxFCr8TNSnehFI5\nShE2PoxL/yX+njVX9i+9vY1F23/+GWZa8Jq09HJtJBc2kgv7SOF3svdLvU8e/zyEjQ8jKjrK7HDi\nFRQEs2ZBmzaw1bqTjgohHER6/C5wJ/YOQ1YNwcvDi/AW4aT1Smt2SPH6/XdjNs/16yF7drOjEUIk\nhvT4LcrTw5OeL/Tk6q2r1J1a17Lz+rz5JjRrZgzzvH3b7GiEEM4ihd9F0nqmpW+lvhy8dJA3p7/J\n4/7SMrN/+ckn4OcH3bqZFsIDpJdrI7mwkVzYRwq/C/mk9WFg2EA2ntpIi1ktHlv8zeLhAZMmGSt4\nTbDu8sJCCDtIj98El29cps+yPtQpVIfv6nxndjjx2rkTKleGhQuhVCmzoxFCJER6/ClEQPoABlYe\nyKy9s+i6sKvZ4cSreHH49luoXx8uXDA7GiGEI0nhN0mgTyADwwYyYfsE+i3r98BrVulfNmwIjRrB\nW29BTIw5MVglF1YgubCRXNhHCr+JgvyC+CTsE8ZsHMOwVcPMDideQ4YYff+PPjI7EiGEo0iP3wKO\nXD7CgH8G0D+sP53LdjY7nEdcvAjPPw/Dhhl/AQghrEN6/ClUvoB8fFzpY/ov78+Pm340O5xHZMli\nTOvQoQPs2GF2NEIIe0nht4iCWQrS+4XedF3Uld4/9TY7nEeULAlffgmvvw6XL7vuuNLLtZFc2Egu\n7COF30KKZStGjwo9+Gr9V0zYbr1B9E2bQp068PbbKWvBdiHEg6THb0Fbz2xlxJoRfFv7W5qWaGp2\nOA+4fRuqVIGqVY11e4UQ5pIefyoRGhRK9wrdaTevHZO2TzI7nAekSQNTpxrz+C9aZHY0QojkkMJv\nQTvW7zCKf/nutJ3X1nLFP2dOmDwZmjeHCCevMSO9XBvJhY3kwj5S+C0sNId1i39YGHzwgXGTV3S0\n2dEIIZJCevwpwNbTWxmxdgRja4+lSYkmZodzT2ws1KsHefPC6NFmRyOEe5IefyoVmiOUbuW70WZe\nGyZvn2x2OPd4eMD48TBvntH3F0KkDFL4LWjH+kfvkiqZoyTdK3Snzbw2lmr7BATAjBnQsSPs2eP4\n/Usv10ZyYSO5sI8U/hTk7miftvPaMm7LOLPDuadkSWM6hwYNIMqaywoLIe4jPf4UaPvZ7QxfM5zh\nVYfT+vnWZodzzzvvwM2bxkIuKkkdRyFEckmP302UyF6CPhX70HNJT75c96XZ4dwzZgzs3m3M4y+E\nsC4p/BYUX4//YUWzFqX/S/0Z9M8gBq8Y7IKonix9eqPf/8knsH69Y/YpvVwbyYWN5MI+UvhTsAKZ\nC/DJS5/w5bov6busr9nhAFCgAPzwA7z5pqzcJYRVSY8/FTh59ST9wvvRtERTvqjxhdnhANCjB2zf\nbiza7iGnF0I4jfT43VSwfzCDKw9mys4ptJ7TGiv8Mh8yBP77D4YONTsSIcTDpPBbUGJ6/A8L8gti\nyMtDmH9gPs3/bG568ffygilT4JtvwJ52rPRybSQXNpIL+0jhT0UCfQIZUmUIq46vov60+sTcMWmF\n9DjBwTBunDF//9mzpoYihLiP9PhToahbUXyy4hNyZMjB/Lfnkz5NelPj+fhjY5TPggXg6WlqKEKk\nOtLjFwBk8M7Ap5U/5cqtK1T8pSKRNyJNjWfAAGMBl8HWGHUqhNuTwm9ByenxP8zby5u+lfrik8aH\nsj+V5fS10w6ILHnu9vvHjoVly5L2Xunl2kgubCQX9pHCn4p5eXrRvUJ38gfkp8xPZTh46aBpseTI\nAb/9Zqzbe+aMaWEIIZAev9uYsH0Cy48uZ2HThZTKUcq0OPr3h5UrYfFi6fcL4QjS4xcJalaiGa8V\neo0q46sQfiTctDj69TMmcBs40LQQhHB7UvgtyBE9/vi8VuQ1Woa2pO7Uuvyx5w+nHONJPD2N2Tt/\n+sk4638S6eXaSC5sJBf28TI7AOFaL+d7Gd+0vrSc1ZLT107Tvkx7l8cQFAQTJhj9/o0bjcXbhRCu\nIz1+N7X3wl6GrhrKe6XeY1jVYSgTJtAfONAY5bNkiTHyRwiRdNLjF4lWJLAIQ18eysQdE2k8szF3\nYu+4PIY+fSBNGmOcvxDCdaTwW5CzevwPy+mfk8+qfMaWM1uoPL4yUdGuXTfR0xMmTjSmdVi4MP5t\npJdrI7mwkVzYx+7Cr5SqqZTaq5Q6oJTqmcA2o+Ne36aUKmnvMYXjZEyXkU/DPuVmzE3K/FjG5Td6\nZc9uXOxt2RJOnnTpoYVwW3b1+JVSnsA+oCpwEvgXaKy13nPfNrWADlrrWkqpssBXWuty8exLevwm\nio2NZeymsWw/u51FzRbxdLanXXr8wYONs/5ly6TfL0RSmNHjLwMc1Fof1VrfBqYCrz20TV1gPIDW\nej2QSSmV3c7jCgfz8PCgXel2VM5XmYq/VmT5keUuPX7v3sbSjf36ufSwQrgle8+tgoGI+56fAMom\nYptcwCMT9V6+ednOcFKHPf/uoWjpoqYcu3pIddJ7pafu1LoMDBtIk2eauOzYX4xVVK+UheLPX6VK\n9WgAVq9czQuVXnBZDFYmubCRXBjSeKZJ1vvsLfyJ7RM9/GdIvO+bNHASQbmCAPD186VAsQKElgsF\nYOu6rQBu8fyCzwUu7Lpg2vEbF2+MxzEP+o/rz+IXF1O3UF32b94PQOHnCgOwb9M+pzxv3P8l2r5X\njTe7fY5f5hsAnN973mnHS0nPI/ZFcD7recvEY+bzJYuXsOr4KsvE48rn+zbtY+2ctQBkDc5Kctjb\n4y8HDNBa14x73huI1Vp/dt82Y4FwrfXUuOd7gZe01mcf2pe2yj0FwnDg4gFqTKxBUIYgupbrSlqv\ntC457owfCrBheXaG/LYGrzTyPSFEQk5cO0G70u1c3uPfCBRUSuVVSqUFGgEPX6GdDTSHe78oIh8u\n+sKaCmYpyObWm1Eoei3rxeUbrmnF1f/fQXz9bjPxqyIuOZ4Q7sauwq+1jgE6AAuB3cA0rfUepVRr\npVTruG3mA4eVUgeB74F2dsac6llpjHKmdJn4p+U/VMhVgW6Lu3H48mGnH9PDAz78bCsr5gUz/fuI\nJ7/BTbjq/o6UQHJhH7vH8Wut/9ZaF9ZaF9BaD4373Pda6+/v26ZD3OvPaq0323tM4Vpenl6Mf308\nXcp3oe/yvqyNWOv0Y/oHRNNt5Cb++LkA506au3SkEKmNpebqsUosImGz9syi5V8tea3wa7xR7A2n\nH+/Pn/OzelFOhk5YTZq08v0hxP3M6vELN1OvaD3CW4az5MgSRq0bRcydGOce753DZMp8i/Ejizn1\nOEK4Eyn8FmSlHn98QoNC2fL+Fq5FX6PXsl5cunHJacfauWEHHwzbyrolQaxZFOS046QE0te2kVzY\nRwq/SJZsGbKx5p01VMhdga6LurLn/J4nvymZMmS8Tc9RG/luQAlOH/dx2nGEcBfS4xd2G7NhDL2W\n9qLZM814peArTjvO3Al5Wfpnbj6bspq03rFOO44QKYX0+IVp2pdpz/y35zNzz0y+2fCN0+b2r930\nKNlz/8fPw1w7gZwQqY0Ufguyeo8/PpXyVGJL6y2c/+88vZf1dtjNXvf3cpWCjp9uY+uaQP6Z637r\nNUpf20ZyYR8p/MJhcvrnZN3/1vF8jufpurgrey/sdfgxfP1i6DlqEz8OLs6Jwxkcvn8h3IH0+IVT\njF4/mj7L+tDkmSbUKVjH4ftf+Hse5k7Ix+e/r8I7veuXjRTCCqTHLyylU9lOLHh7AbP3zebzNZ8T\nHRPt0P1Xb3icfEWvMnZQcYfuVwh3IIXfglJijz8+L+R5gR1tjF7sh4s+JOJq0ufdSaiXqxS07b+d\nfVsDWPpHLrviTCmkr20jubCPFH7hVIG+gYS3DKfJM03ouaQn4UfDHbbv9L536PXVRn4dUYxj+/0c\ntl8hUjvp8QuXmX9gPs3+bEaZnGV4v9T7eHk6ZnHdpX/mYsYPBflixgrS+0q/X7gP6fELy6tVsBbb\nWm/j3H/n6LGkB+eizjlkv1VeP0Gx5y7ybf8SyLmDEE8mhd+CUkuPPz65MuZi/bvrqRFSg66Lu7L2\nxOOneE5sL/f9j3dy7IAfC6c95YgwLUn62jaSC/tI4Rcu5+XpxTe1v+HHV3/ku43f8e3Gb7l957Zd\n+/ROF0vPUZuY+FVhDu3K6KBIhUidpMcvTBVxJYIGvzfg/PXzdK/QndwZc9u1v5XzczJhVBG+nLkC\nXz/nThkthNmkxy9SpNwZc7P23bW8Vfwtei7tyfwD8+3aX6VapyhV8RyjPwqVfr8QCZDCb0Gpuccf\nH08PT4ZWHcqct+bw176/GLpyKFG3ooDk9XLf7bWb86fSM2dCPkeHairpa9tILuwjhV9YRli+MHa3\n301mn8x0XtiZned2Jms/adLG0mPUJqaPLciezQEOjlKIlE96/MKSxmwYQ++lvalZoCZvP/M2Hirp\n5yj/Ls/Gd5+UYOSMFQQEOnbKCCGsQHr8IlVpX6Y96/63jl3nd9FtcTciriR9uofSlc9RpX4EIz58\njjsxSfq5ECJVk8JvQe7W409IsazFGF14NPUK16Pn0p78uedPkvpX4Vvt95EmbSy/fVHESVG6jvS1\nbSQX9pHCLyzN09OTEdVHsKjpIsKPhdNneR/OXU/8Hb+entD1882sXpCTNQtzODFSIVIO6fGLFOPG\n7Rt0XtCTkqxkAAAYOUlEQVSZabum0bxEc6qHVE/0ew/syMgn75dl2KQ15Mof5cQohXAd6fGLVC99\nmvT88OoPTG84nZl7ZjJoxSAib0Ym6r0Fn7lC8657GNrxef6L8nRypEJYmxR+C5Iev018uageUp19\nHfYRkjmETgs6JXqq5+pvRFCk5CW+/jhl3twlfW0byYV9pPCLFMnP248pDabw62u/MmnHJAb+M5CL\n/1184vta993J2Qgf/hqX3wVRCmFN0uMXKd61W9doP789f+37i8ZPN6ZWwVoolXDL8+yJ9HRrVIme\nozZSvPQlF0YqhGNJj1+4LT9vP357/Tf+ePMPFhxaQO9lvTl17VSC22fPdYMPh23h8y7PcfGstwsj\nFcIapPBbkPT4bZKSiyr5q7Cvwz6q5KtC10VdmbZrGrE6Nt5tS1U6T83GR/nsg+e5HZ0ybu6SvraN\n5MI+UvhFquLt5c3oV0azvMVyNp/eTJdFXThw8UC8277Z5gAZ/G/z6/BiLo5SCHNJj1+kWndi7zDo\nn0GMXDeSSnkq0fLZlqRPk/6BbaKupKFLw0q81W4/L9c7YVKkQiSP9PiFeIinhycDKg9gW5ttRN+J\npt38do8M/cyQ8TYfj9nAL58VY//2TOYEKoSLSeG3IOnx2zgiF/kD8rOk+RK+q/0dU3ZO4ePlH3Py\n6sl7r+cpGEWHQdsZ1ul5Lp2z7sVe6WvbSC7sI4VfuI2GTzfkYMeDvPTUS3Rf3J1xW8dxK+YWAOWq\nnqHaG8cZ1ul5bkfLj4VI3aTHL9zSzrM7eW/uexy+fJhmzzTjpbwvERsLwzo/j1/GaDoM2s5jbgUQ\nwhKkxy9EEhTPXpy1767li+pfMHnnZHot7cXRyMN8MHQL+7YFMH9yXrNDFMJppPBbkPT4bZydi7dL\nvM3hzoepW6guH4d/zE+7vqTTl0uY+m0hdqzP4tRjJ5X0tW0kF/bxMjsAIcyWzisdQ6sOpW3ptnT6\nuxMDd7Si7EctGdytEnnLD0N53yStSkudmnUoXaG02eEKYTfp8QvxkBXHVtB8VHOObzyJrhpz7/NB\n64N4r+F7UvyFZUiPXwgHefGpFykcWfiBog9wpuwZ5i6Ya1JUQjiOFH4Lkh6/jVm5uKVvxfv5a7ev\nuTgSG+lr20gu7COFX4h4eKv4b+Q6fOkwo9aP4sL1Cy6OSAjHkcJvQWFhYWaHYBlm5aJTk06EbAl5\n4HMeM0Po8+pPZPXJSscFHflu43dE3kjc0o+O8EzZZ1x2LKuTXNhHLu4KkYB5i+fx9ZSvuRl7k3Qe\n6SidryO/jK3N2rVww2cfPZf0ZMnhJVTJV4VGTzfCz9vP7JCFm0nuxV0p/BYUHh4uZ/1xrJaLzz+H\nCRNg1Srw84NtZ7bRc0lPVkespkZIDd4o+ga+aX2dcuwd63fImW4cyYXB5aN6lFKZlVKLlVL7lVKL\nlFLxTm2olDqqlNqulNqilNqQ3OMJYQVdu0Lp0vD223DnDjwb9CwLmi5gWfNlXLxxkdbzWjNu6ziu\n3TLvIrAQT5LsM36l1HDggtZ6uFKqJxCgte4Vz3ZHgOe01o9d3FTO+EVKER0NNWtCyZIwcuSDr607\nsY5+y/ux9sRaKuetTMOiDcmUXqZ7Fs5hxjj+usD4uI/HA/Ues61MdyVSjbRpYcYMmDsXfvjhwdfK\n5SrHomaLWNlqJRpN2/lt+fbfb2UUkLAUewp/dq312biPzwLZE9hOA0uUUhuVUu/ZcTy3IeP4baya\ni8yZjcLfrx8sWvTo66FBocxpPIeN72/Ez9uPDgs6MHLtSI5FHkv2MWXsuo3kwj6PnatHKbUYCIrn\npT73P9Faa6VUQn2aF7TWp5VSWYHFSqm9WuuV8W3YsmVL8ubNC0CmTJkIDQ29d2HvbgGQ5+71/C6r\nxPPw8xkzwqhfH4YMCadAgfi3/73h70ydM5UJ2yfQ+3RvCmYuSOiNUApkLnDvAuXdQva454f3HE7S\n9qn5+eE9hy0Vjyuf71i/g6V/LgUgXdZ0JIc9Pf69QJjW+oxSKgewXGtd5Anv6Q9Eaa1HxvOa9PhF\nijRzJnTuDKtXw1NPPX7byJuRfLHmC8ZuGouftx+1C9bm5bwv4+nh6ZpgRapiRo9/NtAi7uMWwKyH\nN1BK+Sil/OI+9gWqA/I3mkhVGjSA7t2NC76XHjuEATKly8TAlwdysstJelTowYKDC3h/7vtM2jGJ\nqzevuiZg4fbsKfzDgGpKqf3Ay3HPUUrlVErNi9smCFiplNoKrAfmaq3j6YiK+1m1r22GlJKLzp2h\ndm147TW4efPJ26fxTEPr51uzr8M+xtUbx4X/LvD+3PcZsWYE+y7si/c90te2kVzYJ9nz8ccNz6wa\nz+dPAbXjPj4MhCY7OiFSkOHDoUkTaNYMpk0Dj0ScVimlqFWwFrUK1uLw5cOMXDOSIauGEOgTSI2Q\nGlTOW5k0nmmcH7xwK3LnrhAOdOsW1KhhjPH/8stk7iPmFr9u+ZUxG8cQcSWC8rnLU6tALfJmyuvQ\nWEXKJ1M2CGERkZFQsSK88w506WLfvv49+S+j14/mr31/EZQhiMp5K1MlXxXSp0nvmGBFiiYLsaQi\nKaWv7QopMReZMsHff8OoUTB+/JO3f5zSwaWZUH8CZ7qeoU7aOmw+vZl3Zr/DyLUj2XF2B+56siQ9\nfvvImrtCOEHu3MaNXZUrg78/vP66ffvzSetD3cJ1+SLsC/Ze2Ms367/hqw1foVCUy1WO6iHVye2f\n2zHBi1RPWj1CONHmzcYwz8mToeojQyHso7Vm0aFF/LzlZxYcXEBW36xUyFWBavmrEZA+wLEHE5Yk\nPX4hLGrFCnjjDZg9G8qVc84xbsXcYvru6YzbOo41EWvIF5CPMjnLEJY3jMzpMzvnoMJ00uNPRVJi\nX9tZUkMuXnwRxo0zxvjvsKM1/bhceHt507REU5Y0X8LprqfpWLojRyOP0nZeW7ov7s7vu37n3PVz\nyT+4xUiP3z7S4xfCBWrVgq++Mto+//wDBQo471gZ02WkTek2tCndhqu3rjJ913Sm7ZpGpwWdyJEh\nB6FBoVTIVYECmQuglEyc646k1SOEC/3wAwwdCitXQq5crj32jds3+HPvn8zaO4tlR5ah0TyT7RnK\n5CxD6eDSpPNK3oRfwjzS4xcihRg+HH75BZYvhxw5zIlBa826E+uYuWcmfx/8m8OXD1MgoABPZ32a\n53M+T6EsheSvAYv7d82/zJg7gz1/7JHCnxpYbZ1ZM6XWXAwaZIz0CQ+H7AmtZPEQZ+bi9LXT/LXv\nLxYeWsiaiDXcuH2DwlkKUzxbcUrnLE3ujNYaKurua+7+u+Zffpz+I2fKnoEBJLnwS49fCBP07Qux\nsfDyy8aZf7Zs5saTwy8HbZ5vQ5vn2wCw69wu5u2fx+LDi/lj7x94eXgREhBCoSyFKJGtBIWyFJKp\npE005+85RtFPJjnjF8JE/frBn3/CsmWQNavZ0cRPa82WM1tYengp/xz7h82nNxN5M5J8AfkoEFCA\nwlkKUyxrMbL4ZDE71FQr8kYkO8/vZO+FvRyJPMKu33cRGxZrvDgg6Wf8UviFMJHWxtn/7NlG8Q8M\nNDuixIm4EsHSI0tZcXQFm89sZv/F/aRPk548/nl4KuNTFMpSiKKBRQn0TSFfkIVc+O8C+y7u49Cl\nQxy7coyjkUe5Fn2NgpkL8lyO5yifuwIDOv7G6RqrjDcMkMKfKqTWvnZyuEMutIY+fWDePFi8OOG2\nj5VzcSf2DjvO7WBNxBo2nNzAtjPbOHDpAB7KgxwZcpDDLwe5/XPf+yvB3juLU0OPP/JGJEeuHOFY\n5DGOXTnGyWsnOXn1JHf0HUICQng669OUylmKF/O8yHM5nsPL0wutjbUfFiybR0y+zhx5/pD0+IVI\niZSCwYMhbVrjZq+lSyE42OyoksbTw5PQoFBCg0JpV7odYLSIDl8+zObTm9l2dhs7zu7g912/c/zK\ncTw9PMnqk5VAn0Cy+mQlR4YcBPsHk8c/D4G+gXiolH9v6Z3YO5z/7zynrp3izPUznLp6irPXz3Lu\n+jnOXT9HrI4l2C+Y/AH5KRlUklYlW1E2uCwhASHxjqi6cwfat4etW+Hf1bVZtQG+nvI1C1mY5Njk\njF8ICxkxAr77zij++fKZHY1zaK05duUYu8/tZt/Ffey/uJ9Dlw9x/MpxTl07xY2YG/h7+xOQLoCM\n3hnJmC4jmdNnJkv6LGRKl4mM6TKSKV0mAtIF4JvG1+XDTm/G3OTyjctcunGJSzcvcfnGZSJvRnL1\n1lUib0Zy+abx2pVbV/BN40tW36zk9MtJgYACFAksQtGsRXkm2zPkyZgn0bHfvg0tW8LJkzBnDvj5\n2V5TSkmrR4iUbswY+Owzo+1TuLDZ0bje9ejrHL58mKNXjnL8ynFOXj3JiasnOHXtFJduXLpXZKOi\no7gdexvfNL74pPHB29ObtF5pjX890957eHl44aE8UErhQdy/ygMPPLij7xCjY4jVscTciftXxxAT\nG8OtmFvcjLn5yCNWx+Lv7U/GdBkJSBdAoE8g2Xyzkd03O8H+weTLlI+QzCGEBITgm9bX7nzcvAlv\nvWUU/xkzIP1DSzFI4U8lrNzLdTV3zcWvvxp9/wULoEQJ43Pumov43M3Ff9H/cfb6WS78d4Frt65x\n7fY1om5FERUdxfXb17kefZ3oO9HExMZwR98hVsdyJ/bOvY89PTxJ65GWNJ5pHvhl4e3pbRT3uL84\nMqXLRKZ0mcjikwW/tH4u+yvj+nWoVw8CAmDiRKMd+LDkFH7p8QthQa1agY8PVKtmnOVVqmR2RNbk\nk9aHfGnzkS8g9fXFIiOhdm3jr74ffwRPB942IWf8QljYokXw9tvGHD/2LuYiUo6ICGNiv6pVYeRI\n8HjMte7knPGn/EvnQqRi1asb7Z727WHsWLOjEa6wYwe88AK0aAFffPH4op9cUvgtKDXMQe8okgt4\n7jljNs9Bg8Lp188Y9+/uUuv3xfLlUKUKDBsG3boZQ32dQQq/EClASAh88w3Mnw/vvQfR0WZHJBxt\n6lRo1Mj4t0kT5x5LevxCpCBRUdC4sfHvzJmQWVZVTPG0Nobvjhlj3L19dxRXYkmPX4hULkMGmDUL\nnn8eypaFffvMjkjY4+ZNaNYMpk+HtWuTXvSTSwq/BaXW/mVySC5s7ubC09O4w7d3b2OKh8WLzY3L\nDKnh++L0aXjpJePGLFevyCaFX4gU6p134PffjTPGMWPkom9KsnEjlCkDr75q9PR9fFx7fOnxC5HC\nHTpk3N1ZqpQxz4+ri4hImkmT4IMPHHdvhvT4hXBDISGwbp2xole5cnDggNkRifjcvAlt2sCAAcYk\nfGbekCeF34JSQ//SUSQXNo/Lha8v/PYbtG1r3Pwza5br4jJDSvu+OHzY+H+5eBE2bXLdRdyESOEX\nIpVQyij8c+carYQePWS8vxXMmmX8JdaypXFNxt/f7Iikxy9EqnThgjHR2+nTRk/ZHad3NtuNG8bI\nq1mzYNo0Y/itM0iPXwgBGGv3zp4N//sfVKxozPMj51Wus3mzMdXG6dPGx84q+sklhd+CUlr/0pkk\nFzZJzYVSxsXElSuNaX1few3OnXNObK5m1e+LmBhjGc2aNeHjj42hmla8u1oKvxCpXJEixl2hTz8N\nzz5rFCM5+3e8AweMG+qWLzcu4DZp4rxJ1uwlPX4h3MjatcYkb3nzwrffQp48ZkeU8kVHG3dSf/kl\n9O0LHTs6ZyrlhEiPXwjxWOXLGz3ncuWMG75Gj4Y7d8yOKuVauRJCQ41fqJs2QefOri36yZUCQnQ/\nVu1fmkFyYeOoXKRNa/SfV60ylnWsUAHWr3fIrl3G7O+LixeNC+eNG8OgQTBnDjz1lKkhJYkUfiHc\nVJEiEB5ujP2vX99Y4jEiwuyorO3WLWNVrKJFIX162L0bGjSwbi8/IdLjF0IQFQXDhxuTvbVrBz17\nGlNAC4PWxlj8jz4yLpIPG2b8awXJ6fFL4RdC3BMRYRS3ZcugVy/jQnC6dGZHZa7wcOMu6NhY4yJu\n5cpmR/QgubibSpjdv7QSyYWNK3KROzdMmGBM+7BkCRQoYCz5ePOm0w+dJM7OhdbG1//ii0Yvv3Nn\n2LDBekU/uaTwCyEeUbIk/PWXMd3AggWQL59xY9LFi2ZH5lx37hhLWlaoYAzLfP992LvXuP6REkbr\nJJa0eoQQT7RzJ4wcafwiaNLEuCBcvLjZUTlOZCT88gt8/TXkzGlMcle/vrHamdVJq0cI4RTFi8Ov\nv8KuXZAlC9SoYUwz/NtvcP262dElj9ZG/75ZM+OGto0bjQu4q1dDw4Ypo+gnV7ILv1KqoVJql1Lq\njlKq1GO2q6mU2quUOqCU6pnc47kT6WvbSC5srJCLnDlh4EA4dgy6dzcKZXCw0QqZN89100AnNxda\nw5YtxqyZBQpAhw7GZGoHD8LkycZyiO7AnjP+HcDrwIqENlBKeQLfADWBYkBjpVRRO47pFrZu3Wp2\nCJYhubCxUi68vIzlHufNg/37jZ744MGQPTs0agQTJzr3ekBScnHrlrEgfZcuUKiQMe5ea5g+HXbs\nMNo6gYHOi9WKvJL7Rq31XjD6S49RBjiotT4at+1U4DVgT3KP6w4iIyPNDsEyJBc2Vs1FtmzQvr3x\nOH0a5s837ghu1w7y54eXXjIeZcsafzE44manx+Xi+nVjBM6qVUbbZs0ao1X1yivGBHWlSqW8G64c\nLdmFP5GCgfvvBTwBWGxmaiGEo+TIAe++azxu3zbmrwkPh59+gtatjZExpUoZhbhgQaPdEhICQUHg\n7Z3442gNZ8/C0aPG48AB2L7deBw/bsyf88ILxrTUkyYZ1yWEzWMLv1JqMRAUz0sfaa3nJGL/Mkwn\nGY4ePWp2CJYhubBJablIk8aYDK5cOeNmMK3h5Enjl8Hu3cb8QJMmwaFDxjoBPj7GXw/+/sZ0COnS\nGfu4fdt4REfD5ctGC+n8+aN8+61xUTZvXuOXR7160K+fsdpYmjRmf/XWZvdwTqXUcqCr1npzPK+V\nAwZorWvGPe8NxGqtP4tnW/klIYQQyZDU4ZyOavUkdNCNQEGlVF7gFNAIaBzfhkkNXAghRPLYM5zz\ndaVUBFAOmKeU+jvu8zmVUvMAtNYxQAdgIbAbmKa1lgu7QghhIsvcuSuEEMI1XHrnbmJu5lJKjY57\nfZtSqqQr43OlJ+VCKfV2XA62K6VWK6VKmBGnKyT2Jj+lVGmlVIxSqr4r43OlRP6MhCmltiildiql\nwl0cossk4mckUCm1QCm1NS4XLU0I0+mUUr8opc4qpXY8Zpuk1U2ttUsegCdwEMgLpAG2AkUf2qYW\nMD/u47LAOlfF58pHInNRHsgY93FNd87FfdstA+YCDcyO28Tvi0zALiBX3PNAs+M2MRcDgKF38wBc\nBLzMjt0JuagElAR2JPB6kuumK8/4793MpbW+Ddy9met+dYHxAFrr9UAmpVR2F8boKk/MhdZ6rdb6\nStzT9UAuF8foKon5vgDoCMwAzrsyOBdLTC6aADO11icAtNYXXByjqyQmF6cB/7iP/YGL2riumKpo\nrVcClx+zSZLrpisLf3w3cwUnYpvUWPASk4v7vQvMd2pE5nliLpRSwRg/9N/FfSq1XphKzPdFQSCz\nUmq5UmqjUqqZy6JzrcTk4kfgaaXUKWAb0NlFsVlNkuums+/cvV9if1gfHtaZGn/IE/01KaUqA+8A\nLzgvHFMlJhejgF5aa62MOUJS69DfxOQiDVAKqAL4AGuVUuu01gecGpnrJSYXHwFbtdZhSqkQYLFS\n6lmt9TUnx2ZFSaqbriz8J4Hc9z3PjfGb6XHb5Ir7XGqTmFwQd0H3R6Cm1vpxf+qlZInJxXPA1Lh5\noQKBV5RSt7XWs10TosskJhcRwAWt9Q3ghlJqBfAskNoKf2JyUQEYDKC1PqSUOgIUxrh/yJ0kuW66\nstVz72YupVRajJu5Hv7BnQ00h3t3/UZqrc+6MEZXeWIulFJ5gD+AplrrgybE6CpPzIXWOr/WOp/W\nOh9Gn79tKiz6kLifkb+AikopT6WUD8bFvN0ujtMVEpOLvUBVgLiedmHgsEujtIYk102XnfFrrWOU\nUndv5vIEftZa71FKtY57/Xut9XylVC2l1EHgOtDKVfG5UmJyAfQDAoDv4s50b2utU91s4YnMhVtI\n5M/IXqXUAmA7EAv8qLVOdYU/kd8XQ4BflVLbME5ie2itL5kWtJMopaYALwGBcTfN9sdo+SW7bsoN\nXEII4WZk6UUhhHAzUviFEMLNSOEXQgg3I4VfCCHcjBR+IYRwM1L4hRDCzUjhF0IINyOFXwgh3Mz/\nAZRHjy2qisvqAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# fix nodes\n", "h = 1\n", "x = np.linspace(0, h, 3)\n", "\n", "# find interpolant\n", "V = np.array([\n", " 1+0*x,\n", " x,\n", " x**2\n", "]).T\n", "a, b, c = la.solve(V, f(x))\n", "interpolant = a + b*plot_x + c*plot_x**2\n", "\n", "# plot\n", "pt.plot(plot_x, f(plot_x), label=\"f\")\n", "pt.plot(plot_x, interpolant, label=\"Interpolant\")\n", "pt.fill_between(plot_x, 0*plot_x, interpolant, alpha=0.3, color=\"green\")\n", "pt.plot(x, f(x), \"og\")\n", "pt.grid()\n", "pt.legend(loc=\"best\")" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.161228524153\n" ] } ], "source": [ "true_val = int_f(h)\n", "error = (f(x).dot(simpson_weights) * h - true_val)/true_val\n", "print(error)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Compare the error for $h=1,0.5,0.25$. What order of accuracy do you observe?" ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 0 }