{ "cells": [ { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "---\n", "title: Numerical Integration\n", "date: 2019-09-01\n", "\n", "# Put any other Academic metadata here...\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the integration of the function $f(x) = \\exp \\left( x^2 \\right)$ between $[0,1]$.\n", "\n", "It can easily be shown that\n", "\n", "\\begin{align}\n", "\\int_0^1 e^{-x^2} \\, \\mathrm{d}x = \\dfrac{\\sqrt{\\pi}}{2}. \\nonumber\n", "\\end{align}\n", "\n", "First import the tools needed" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.integrate as spi\n", "from IPython.display import display, HTML\n", "# some formatting for the graphs \n", "\n", "display(HTML(\"\"\"\n", "\n", "\"\"\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next define the range to plot the function over and a number of point to evaluate the function." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "x0 = -0.5\n", "x1 = 1.5\n", "n = 100" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So we have a range and a function. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "x = np.linspace(x0, x1, n)\n", "y = np.exp(-x**2)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "f = lambda x : np.exp(-x**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let $m$ be vector whose elements are the number of evaluations of the integral. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "m = np.array([4,8,16], dtype=int)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute and plot the successive approximations for the integral, as well as keeping track of the errors.\n", "\n", "We will use the function [trapezoid](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.trapezoid.html#scipy.integrate.trapezoid) to compute the integral according to\n", "\n", "\\begin{align}\n", "\\int_a^b f(x) \\,\\mathrm{d}x \\approx \\sum\\limits_{k=1}^{M} \\dfrac{f\\left( x_{k-1}\\right) + f\\left(x_{k}\\right)}{2} \\Delta x_{k} \\quad \\mbox{where} \\quad x_0 = a \\quad \\mbox{and} \\quad x_M = b. \\nonumber\n", "\\end{align}" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAFACAYAAADku5hoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABzgklEQVR4nO2dd3hcxdW436PeLctFstzkIndcZRtsY2RjMBgCJiRgSoAEAiQh9fuSwAchBEJCQvKDJCQBQjEBginBxBBTDTK4W+4Y9967ZEuy+vz+OHfR9XpXu7K1Wq007/PMs/fOzJ17Zu7sPXfaGTHGYLFYLBZLpBAVbgEsFovFYmkIVnFZLBaLJaKwistisVgsEYVVXBaLxWKJKKzislgsFktEYRWXxWKxWCIKq7gsFovFElFYxWWxWCyWiMIqrlaCiKwVkXw/YdNF5Nchuu92EZkUirQtFkvrxCquIBCREperFZGTrvMbwi1fMBhjBhpjCs7kWkf5ePK831F0KY0sYjBytBURIyILvfyfFJHHQnC/XBEpF5GXzvD6AhE5JiLxjS1bUyAiGSIyU0RKRWSHiFxfT9wSL1cjIn9xhb8kIvtE5LiIbBSR21xh/UXkYxEpFpHNInKVV9qBwnNEZLZT1vtF5AkRiXGFTxORdU4+tojI+T7kP+1ZB5GnAucaT/gGV1i8iDzrlNsJEVkpIpcGm3YwcvsLD6I87hKRQhGpEJHpPsqi3vJuFhhjrGuAA7YDk+oJjwm3jGeQp+nAr4PJM5AFrAIebozyaqCcE4B9QDGQ5fJfDNwcgnL5APgMeOkMrs0BaoCjwNdD9NxCWteAV4BXgRRgnFPuA4O4LgUoAca7/AYC8c5xP2A/MAKIATYCPwGigYlAKdDHk8f6wp04s506nODUzzXAD5ywi4AdwLnoh3pnoHNDn7WfPBUAt/mJnww84NSDKOBy4ASQE2Ta9cpdX3h95eGEfxWYCvwdmO5dpwKVd3NwtsXVCDgtkp+LyGqgVERiRORu5yvohIh84f5qceLf4/gfE5HnRSTBFZ4tIv8WkUMisk1EfuD4X+v1lVYhIgVOWH/nC7BItFvwCh8yTnKOh4nIcke2V9EKHhTGmP3A+8BQV9pGRHq7zv12PfrLW5AMBQqBD4ErnfSigXOAFQ1IJyAiMg0oAuacYRI3AYvQF8jNXmkHev5+w/3UNZ/PXkR6ichRERnunGc75Z4fRP6TgauBXxhjSowx84BZwDeCyPvVwEFUEQBgjFlrjKnwnDquF6rEsoHHjDE1xpiPgfmu+wQKB+gBvGaMKXfq53uoogT4FfCgMWaRMabWGLPHGLPHK6/BPOvT8lQfxphSY8wDxpjtzn3fAbahyjqYtAPJXV94feWBMeZNY8xbwBEfsgRT3mHHKq7G4zrgMiDdGFMNbAHOB9qglewlEenkin8DMBn98/YB7gMQkSjgbbRV0xm4EPiRiEw2xrxqjEkxxqSglWsr8IqIxDrXfAB0BL4PvCwifb2FFJE44C3gRSADeB394wSFiHQBLgU2B3uN61q/eQsyiWHASlT+qY5fP7Qer/Nzz3ecF7ov946fa9KAB9GvzjPlJuBlx00WkUyvcJ/PP8jwL+saIPh59saYLcDP0bqXBDwPvGBcXcYi8jcR+ZsP+fsA1caYjS6/VbhegPVwM/BP43zCe92rDFiPtpxn+7legEH1pO8d/jgwTUSSRKQzWj/fcz5q8oAOTpfXbqfbLNElU7DP2meegN+KyGERmV/fB4Hz/PsAawOlHUjuIPLlszwC5K8+Aj2PpifcTb5Ic/jo+nL8vhXgupXAla74d7rCpgBbnOPRwE6va+8BnnedRwHvAH93zs9Hu16iXHFeAR7wlhsYD+wFxBW2gMBdhSVoV4dBv0zTXeEG6O06n+5Jj1O7GQPmLUAZrkaVbFtHnlT0Bb+skZ/xn4CfO8cP0MCuQrRbrQpo75yvB37sVZ4+n38Q9eOUuhbks5+FdhetxumuCyIP5wP7vfy+DRQEuK472kXaw094tFM+9wGxjtsK/Mw5vhioBN534tcb7sTpDywDqp26OB192WY754VAJ6A92np42HVtwGftL09OfU4F4lHlcwLo5eP6WOAj4Klg0g4kdxDhPsvDx71/zeldhQHLuzk42+JqPHa5T0TkJtEB2SIRKUK/WNr7ib8DrYygFTnb3TIA/g9wf7E/jP5hPN1s2cAuY0ytV5qdfciZDewxTi11xQ3EVGNMKpCPtnLa1x/dJ8HkzSeiExz6AyuNMceAJeiXpKcV1iiIyFBUwZ/NZI+bgQ+MMYed83/h1V2I/+cfTLg7LJhn/w+0/v3F1HXXBaIESPPyS0NfzvXxDWCeMWabr0Cj3U/zgC7Ad4wxVWjr+TJUAf8P8Bqw24lfb7jTin8PeBMdV2qPftj8Djjp3PYvxph9zvP4f+iHQEOetc88GWMWG2NOGGMqjDEvoMpjijuOI9+L6Mv/riDTrlfu+sIDlEdAApV3cyEmcBRLkHypCESkO/qyuBBYaIypEZGV6Fegh66u425oKwj0pbTNGJPr6yZOf/x1wEinkuFc21VEolwvsG7oIKs3+4DOIiIu5dUN7doMnElj5orORPoDdd11ZUCSK1oWvit6vXkLwCDnPlud87ec+2cCM/1dJCLvoq0HX3xmjLnUyy8fHVDfKSKgA+fRIjLAGDM8kJBOd801zjX7He94IF1EhhhjVjl+/p4/QYS7Pzrqffaisz8fB54FHhCRfxtjjgbKh3N9jIjkGmM2OX5D8N3V5eYm4JEg0o9Bu0ExxqwGLvAEiMgC4AXPeYDwDDS/TzhKuUJEnkdb/D8Tkd2cWl7u43yCe9bB5sng+o+LJvosWkenuP6vbk5L2xhzrD65A4T7LQ+0FRU4EwGeR7Mg3E2+SHP47yqc5DofAJQDfdGukW+izfbbXPHXoF+dGcA84DdOWDSwHB2bSHTOBwEj0dbFIWCo1/3j0Bf63WjzPh/9Mu7nLaMTdyfwQyfuV9FuraBmFTrnHdCZRkOc8/nony8auAT9IvTVVeg3b074dLy6Llz3vA1VNJ7zHHRA/SgwrhGfbxKqeD3uD8AbQAdXnPrkvM6RqZtXOp8Cfwz0/IOoH97Pot5nj744X3WOn0YH7YMtixlot2MyMJYAswqBMU69SPXy7whMw1EM6NhdKXCFEz4YnSCUBPwvOokh3nV9oHBP/mPQcb+ZwL+csAeBpY4MbdEJEA814Fn7y1O6k48E5743cPpsxyfRCTopDSmvQHIHkS+/5eGExzhy/xZtDSbgmqEaqLybgwu7AJHmCEJxOX4Poy8wTzN+LqcqrnuAL9CX7wtAkuvabPSFsR845lT+SWgffDXajeNx7zrXDHTuUeyke5U/GdGB3RXoC+5VxwWtuBy/vwP/dqW31knvRUf20xRXfXlzwuYA3/YjwxNo14jbbyVQ6+uP34jP+wG8xj0CyPkejoLy8r/GyXNMEM/fb7ifZ+Hz2aMzL/cAGc55Cjqp5gbXtU8CT/rJSwbasi1FP3au9wp/F/g/1/lTwIs+0ungyFcEHEeV8rdd4Y86daHESbO31/WBwoeiU9OPof+314BMJywW+Jtz7/3An4GEBjzr+vK0FK3zRU49vsgV3h1tBZVz6v/1hkBpByN3feH1lYcrn8bLPRBseTcHJ46gliZERLajSuyjcMvSXHBmO64CBhvfXSrNgsaQM9Dzt/XDYqkfO8ZlaRYYYyrRyRfNmkiR02JpydhZhRaLxWKJKGxXocVisVgiCtvislgsFktEEbIxLhF5DjUsedAYc5q5EGeNw5/QRXVlwC3GmOWB0m3fvr3Jyck5K9lKS0tJTk4+qzSaikiRtaisiv3Hy0mPraGoKpqstATSk2LDLVa9RErZgpU1VFhZQ0Njybps2bLDxpgO3v6hnJwxHZ3C/E8/4ZcCuY4bjU6vHh0o0ZycHAoLC89KsIKCAvLz888qjaYiEmR9a8Ue7nlzDe2ravifc6r545oYYmOjue+r5zB1mC/jHc2DSChbD1bW0GBlDQ2NJauI+LTqEzLFZYz5VERy6olyJXWGJReJSLqIdDLG7AuVTJbgqK2F0lIoKan7PXkSKiqgvFxddbW6qip4cFYxx8o6gYFFR2sp2SOUCPxsUxFVV3cmLo4vXWKiuqQkSE6G1FRISYHo6HDn2mKxRAohnZzhKK53/HQVvgM8YtRuGSIyBzV2eVpzSkRuB24HyMzMHDFjxoyzkqukpISUlCbfB/GMOFtZa2vhxIlYjh6Nc1wsRUVxFBXFUlwcy/Hj6k6ciKGkJIbS0hjKyqIxRgIn3ogkJNSQklL9pUtLqyItrYo2bdS1bVtF27aVtG1bSbt2laSnV561smtN9aApsbKGhtYo64QJE5YZY/K8/SNiHZcx5mnUXA15eXnmbJugLanJXVIC27bBjh2wfbv+7toFu3er27tXW0XexMRA+/bqMjKgRw9IT4c2bSAtTVtCqanaKkpOrmspxceri4vTNGJi4NqnF7D/eDkicFvfav6xPgaM0CYmkR8OPpeqqrrWWXW1KtOaGj2uqIDKSigri6akJJoTJ+IpLoaiIti5Ew4d0jjeREdDVhZkZ0PXruq6dYPu3SEnR/PTti2ID/27dNZTdF3+KOv73EW/wifYNfynjLzijjN7QE1ES6qzzQkra2gItazhVFx7ONWQaBfHz+LFiROwYQNs3AibNqnbvBm2btUXu5v4eOjSRV/k558PnTvryz0rCzp1go4d1aWn+36pnwm/uK4797y5hpNVNbTrUE3s/hjio6P59og+5HU/NW5Nzamuurru2J88xtR1WZ44oUrt6FE4eBD27YO1a+GDDzTcTVoa9OoFvXvrb58+UHPwLS4+9BiZyYdYL5DFIdosu4+l0OyVl8ViUcKpuGYBd4nIDHRSRnFrH986fFhfwl98ob/r1sHq1edx+HBdHBFVSr17w9Sp0LOnti5ycrS10bEjRDXxIgfPBIxHZm8ATtA+KZEbBvVlfPfTJ2ZERzdsPMsYVWppaacquZqa0+OVlKgiP3xY3YEDqtiWLoWZM/V6NSg/lbYJx8judpJzkzPo334D2VsWkT3iDrKzG0+hWyyW0BDK6fCvoJaq2zsm+H+JGobEGPMkuvvpFNToZxlqQb1VUFWlLaiVK2HVKli9Gtas0Zesh9RU6N8fRow4Rn5+Fn37Qt++qqgSEvwmHTamDuvM2C6dWbiwgKcuy2+0dEXquiQD0bGjKnFfrbrKSti/Hzp+9DU2H+3N+sN9WFh6MbM2XMqzK27SBF7WluigQTBkSJ0bNEgnk1gsluZBKGcVXhcg3ADfC9X9mwsVFaqYli9Xt2wZfP553bhNfDwMGAAXXQTnnKMvyYEDtbtPBAoK1pOfnxXeTEQQ9Sm6zEzI3rmSy/kQgIK+vyJ/wy85XJZBwb4xfJz6Mjt36hjh9OnaPQnagu3XD4YNg+HDIS9Pj1NTmy5fFouljoiYnBEp1NZq996SJeqWLlWl5Zkc0batvvi+/30YOlRd377BtSYsjcOKrj8lfed9JErll35JCaUkjpnAraP0GZaXq9Lau7du4suOHTBnDrz8sl4jos9u5EgYNUrdkCH6IWKxWEKLfWWeBUeOwKJFsGCB/i5dqpMHQMdk8vLgJz/Rl9uIEToGZcdPwkvm0Dv4BBi681EwsLemIyuyf07WqNsAbV0lJanr0EGVUXW1rmMrL9dJIZs364zHrVvhvffgxRc17bg4/TA591w47zx1Xbv6l8VisZwZVnEFiTH6wpo3T938+TpOBTrZYMgQuPFGGD1aXZ8+TT9JwhIcmUPvYGWnOygtLSB16iYCdcTGxNQtD+jQQVtaFRWqzCordSLIxo26HGH9enjySXj8cb22a1cYM0bd+efD4MF2sbXFcrZYxeWHmhqdOPHZZ3Xu4EENy8jQF9HNN+tvXp6udbK0Hjzr2UAnhfTrpy0yjzLburVu2cJnn8Grr2rc1NQ6JTZ+vLbGm+NkG4ulOWMVl0N1NaxYAQUFMHeutqqKizUsJwcmT9aXzbhx+sVtW1MWN+4uRtA1c6NH13UxHjig458bNuhSh/ff13jx8RovP1/duefqQm8P/y34BX/aOpOpbb7Lb567ix/2vIrL8h9q6uxZLM2KVqu4amu1RfXxx/DJJ/pVfPy4hvXtC9deCxdcoF/FXbqEV1ZLZOLdxThwYJ29xyNHVIF9/rkqtF//Gh58UMfJzj0XJk6EhPRneC3xHSrjdWB0X7TwwLaZAFZ5WVo1rU5xffwxPPDAQNau1YF20PGo666DCRNUWWXZ2eeWECCi3YIJCbperFcvuOwyVWSHD+taPo/71a/AmNuQuBtI7rOcD4ZEcbL3AEz39fxp60yruCytmlanuPbuhY0bU5k6Vb9qJ0xQk0gWSziIjq6zB9m9O1x+uS6fOHgQvv/RDyhdP4qSL0bzn5dzgfOITi5mV//FPFkJkyap8rMzVS2tjVanuK6/Hjp3XsSECflhlsRi8U1srNqY7DPkY/aP+ASA66t+wt8K11PyxWhOfj6G73xH43brpuOvF10EF16oE4cslpZOq1NcUVH2C9USGVyfehV/K51JeZTQpu1J0s+bTcdRs/lGzA20k7tZtEitsfzrX/CPf2i9zstTRTZ5sk76iG3em1BbLGdEq1NcFkukMLjXQ3x3C7x84i0A2lXGcHXitYzpfzegpsK+9S218rFsGRQW6szY3/xGJ3ukpWkrzKPIcnLClxeLpTGxistiacYM7vUQPUseYv/+Av44fIXPOMnJOvt1/HidLXvggFpyWbZMF8rP1ImI9O0Ll1yi7oILTp12b7FEElZxWSwtiKgoXUN21VXqKit1m5wlS1SRPfkk/OlPOrMxPx8uvVQVWW6u7UK3RA5WcVksLZi4uDqDzsboJpyLFmm34vLlamsRdDuYKVNUkU2YYLdxsTRvrOKyWFoJIrpDwaWXqquuVrNUCxeqInv2WfjrX9Wax/jxqsimTLGtMUvzwyoui6WVEhOjm5X276+TPIqLYfHium7FDz+EH/9YW2OXXabKLj/ftsYs4cda3LNYLAC0aQMXXwz33QdvvKHT7L/7XTVX9Y9/qPJq107HxP78ZzUg7Kbg2QeZP3oQ5WvXMn/0IAqefTA8GbG0eGyLy2KxnEZ0tJpC69NHW2PHj+vYmGeD1Pffhx/+sM5sVeeolxn/wZu0qanhBJBRXEPF469QAOTfen+Yc2NpaVjFZbFYApKWpq2xiy/WKfebNukOCoWFOlOxsvIGEuSrjEpawsC5wpWVXenGLmKffg2s4rI0MlZxWSyWBhEVpWvC+vaFW2/VXb83fOPbfFYyns9Kx/Pp6935Ox/QPXY745I/5cR7dt2YpXEJ6RiXiFwiIhtEZLOI3O0j/BYROSQiKx13WyjlsVgsjU9qKgzMXsi9mb/hvZ6X8Mz9L3FPx4fpGreT14uv4dJL1Ybi5Mnwl79oa81iORtC1uISkWjgr8BFwG5gqYjMMsZ84RX1VWPMXaGSw2KxhJ69U68h5eVXiK+Gzh2LGdP2Jb7e/l+s/sq3OZDxIxYs0LGxDz7Q+Dk5OsljyhRdN5aSElbxLRFGQMUlIt8HXjLGHGtg2qOAzcaYrU46M4ArAW/FZbFYIpxeF9/PRiD7rdcAOJwcz97Lb2TQFT9iEGozEWDrVt20dckSeP55HR+LjYXzzlMldsklMHiwXTdmqR8xxtQfQeTXwDRgOfAc8L4JdJFe9zXgEmPMbc75N4DR7taViNwC/BY4BGwEfmyM2eUjrduB2wEyMzNHzJgxI6jM+aOkpISUCPnEiyRZq6uhtLSE2NjmL291NdTUlBAf3/xlra2FqqrIkBWgoiI4WauqhM8/b8OyZRksX57Bjh16TUZGBSNGHGPUqKOMGHGMtm2rQiZrJP2/WqOsEyZMWGaMyfP2D6i4AEREgIuBbwJ5wGvAs8aYLfVcE4ziageUGGMqROQO4FpjzMT6ZMnLyzOFhYUBZa6PgoIC8vPzzyqNpiKSZD10CBYuLKBz5/xwixKQAwegtLSAnj3zwy1KQEpKYP/+Anr3zg+3KEGxaVMBubn5Db5u/35tjS1cqFbuT5xQ/8GDdXzskktg7Fi17NFYRNL/qzXKKiI+FVdQY1zGGCMi+4H9QDXQFnhDRD40xvzMz2V7gK6u8y6OnzvdI67TZ4DfByOPxWJpeWRlwde/rq6mBj7/XKfcL10K/+//waOP6szEMWNUiU2eDIMG2W7F1kgwY1w/BG4CDqPK5afGmCoRiQI2Af4U11IgV0R6oAprGnC9V9qdjDH7nNMrgHVnlAuLxdKiiI6GIUPUge45tmgRLFigxoHnzIGf/lStekyYoEps0iTdEdrS8gmmxZUBfNUYs8PtaYypFZHL/V1kjKkWkbuA94Fo4DljzFoReRAoNMbMAn4gIlegrbijwC1nmA+LxdKCSU7WCR6eSR779mlrbMkStan4ms4JoWdPjXPRRWpXsUOHsIlsCSEBFZcx5pf1hNXbQjLGzAZme/nd7zq+B7gnsJgWi8VSR6dOdd2KxsDGjbppZmEhvPSS2lYEGDgQJk7U1tj48ZCerv7Fb7/Nwccep/zqr7LpwYfo+OMf0eYrXwlbfiwNw1rOsFgsEY1InSWPb31LZ4yuWqVdi8uXw1NP6cLnqCjtehzbfQsDVn3I8Gid/VG9dy/7fqHf01Z5RQZWcVkslhZFTAyMGKEOoKJCt2lZvBhWr4an/tOVKvNnhFp6PnKEsaUp5CUWMup3zzHaKq6IwCoui8XSoomP15mIY8bouVw3mjXl57CkbBTzk6by6p5r+eexm2EvDBgI559f5+xkj+aJVVwWi6VVEd8xg1GHlzIqaSlf+b7Q6U9Psab8HFbETmBd99v417+0exGga1cYN07Xj40dC+ecozMeLeHFKi6LxdKqqLrmR8Q9cz9SWQ5AXFQVeRlf8JWHptHmK7qGbPVqnewxbx7MnQuvvKLXpqTAueeqiaoxY2D0aGjbNoyZaaVYxWWxWFoVNWO/QiUQ+9rjAMR0yqbjT+pmFUZHw7Bh6u66S2ct7tiha8jmz9ffhx9WU1wA/fqpMjv3XFVkgwbpOJsldNjitVgsrY6asV9RBbangG7vfa9eM1Iias0+Jweud0wolJToGrKFC3X24jvvwPTpGpaYqBNDRo2CkSPV9expLXw0JlZxWSwWSwNJSdH1YRMdy6rGwLZtqsSWLFH317/qjEbQ/chGjIC8vLoZj927W2V2pljFZbFYLGeJiLaqevasa5VVVam9xaVLdWH00qVqb7G6WsMzMmD48LpuyWHDIDfXTv4IBqu4LBaLJQTExtYppNtvV7/yclizRteVLV+u7k9/gspKDU9MVGv4HjuNgwerS0sLXz6aI1ZxWSwWSxORkFA37uWhshLWr9etXFauVPf66/D003VxuneH7OxB5Ofr5I9Bg9RSSGNu8RJJWMVlsVgsYSQurq5ldfPN6mcM7N6tpqvWrFG3aFHCKV2N0dHQuzcMGFDn+veHPn3UKHFLxioui8ViaWaI6OLnrl3hcmcPjoKCQsaMyWfTJh07+/xz+OILdbNm6fozD9276zR9jw3Hvn1VoXXurDYbIx2ruCwWiyVCiItTi/cDB8K119b5V1aqhfz169WtW6e/8+bpXmYeEhK0lZabe+pvz57QpUvkTAyxistisVginLi4urEvN8bA3r2wYYMqtk2b9HfdOvjvf+smhYBOJsnJqZsd2aOHOs8atnbtms/0fau4LBaLpYUiot2DnTvXrTnzUFOj42ibN8PWrbBli7pt23Qd2rFjp8ZPStIuyO7d1fhwt2513Zldu2qLLSGhafJlFZfFYrG0QqKj6xSRZ2dpN8XFqsS2b1eTV9u2wc6delxYCIcPn35N+/aqwBITB/Hyy9piCwVWcVksFovlNNq0gaFD1fmirExbbLt2qdu9u86tX59AbGzoZLOKy2KxWCwNJilJZyr26XN6WEFBIV265Ifs3i1gYqTFYrFYWhMhVVwicomIbBCRzSJyt4/weBF51QlfLCI5oZTHYrFYLJFPyBSXiEQDfwUuBQYA14nIAK9otwLHjDG9gceA34VKHovFYrG0DELZ4hoFbDbGbDXGVAIzgCu94lwJvOAcvwFcKNJcVgpYLBaLpTkixpjQJCzyNeASY8xtzvk3gNHGmLtccT534ux2zrc4cQ57pXU7cDtAZmbmiBkzZpyVbCUlJaSkpJxVGk1FJMlaU6PyijR/eWtrwZgSoqObv6zGQG1tZMgKUFMTObIaU0JaWkqzWVhbH5H0LmgsWSdMmLDMGJPn7R8RswqNMU8DTwPk5eWZ/Pz8s0qvoKCAs02jqYgkWY1ReceMyQ+3KEExf34BY8fmh1uMoLCyhoYFCwqYMCE/zFIERyS9C0ItaygV1x6gq+u8i+PnK85uEYkB2gBHQiiTJYSIqIuUrRaioqysoSCSZI2ElpbldEKpuJYCuSLSA1VQ04DrveLMAm4GFgJfAz42Afouly1bdlhEdpylbO0BH+u+myWRJCtElrxW1tBgZQ0NrVHW7r48Q6a4jDHVInIX8D4QDTxnjFkrIg8ChcaYWcCzwIsishk4iiq3QOl2OFvZRKTQV79pcySSZIXIktfKGhqsrKHBylpHSMe4jDGzgdlefve7jsuBr4dSBovFYrG0LKzlDIvFYrFEFK1VcT0dbgEaQCTJCpElr5U1NFhZQ4OV1SFk67gsFovFYgkFrbXFZbFYLJYIxSoui8VisUQULVZxiUiGiHwoIpuc37Z+4tWIyErHzXL593As1m92LNjHhVNWERkqIgtFZK2IrBaRa11h00VkmysfQ0Mg4xlb+heRexz/DSIyubFlOwNZfyIiXzjlOEdEurvCfNaHMMp6i4gccsl0myvsZqfObBKRm0Mta5DyPuaSdaOIFLnCmqxsReQ5ETnomJXzFS4i8mcnH6tFZLgrrEnLNQhZb3BkXCMiC0RkiCtsu+O/UkQKm4Gs+SJS7HrO97vC6q07DcIY0yId8Hvgbuf4buB3fuKV+PF/DZjmHD8JfCecsgJ9gFznOBvYB6Q759OBr4VQvmhgC9ATiANWAQO84nwXeNI5nga86hwPcOLHAz2cdKLDLOsEIMk5/o5H1vrqQxhlvQV4wse1GcBW57etc9w23PJ6xf8+un4zHGU7HhgOfO4nfArwLiDAucDiMJZrIFnHeGRAd9tY7ArbDrRvRuWaD7xztnUnkGuxLS5OtTz/AjA12AtFRICJqMX6Bl9/BgSU1Riz0RizyTneCxwEznoxdpCcjaX/K4EZxpgKY8w2YLOTXthkNcZ8Yowpc04XoebIwkEw5eqPycCHxpijxphjwIfAJSGS00ND5b0OeCXEMvnEGPMpatTAH1cC/zTKIiBdRDoRhnINJKsxZoEjC4S3vgZTrv44m7p+Gi1ZcWUaY/Y5x/uBTD/xEkSkUEQWichUx68dUGSMqXbOdwOdQydq0LICICKj0K+WLS7vh53uhMdEpLEtxXUGdrnOfZXHl3GccitGyzGYaxuTht7vVvTL24Ov+hAqgpX1aufZviEiHvufTV2uDbqn0/3aA/jY5d2UZRsIf3kJR7k2BO/6aoAPRGSZ6C4azYHzRGSViLwrIgMdv0Yt14iwDu8PEfkIyPIRdK/7xBhjRMTfvP/uxpg9ItIT+FhE1qAv3UalkWTF+Sp8EbjZGFPreN+DKrw4dP3Ez4EHG0PuloyI3AjkARe4vE+rD8aYLb5TaBLeBl4xxlSIyB1oq3ZiGOUJlmnAG8aYGpdfcyvbiEJEJqCKa5zLe5xTph2BD0VkvdMqChfL0edcIiJTgLeA3Ma+SUS3uIwxk4wxg3y4/wAHnJe852V/0E8ae5zfrUABMAy1UJ8uarEefFu2b3JZRSQN+C9wr9O94Ul7n9PlUQE8T+N3xTXE0j9yqqX/YK5tTIK6n4hMQj8arnDKDfBbH8ImqzHmiEu+Z4ARwV4bAhpyz2l4dRM2cdkGwl9ewlGuARGRwejzv9IY8+UOGq4yPQjMJLTd8AExxhw3xpQ4x7OBWBFpT2OX65kOjjV3BzzKqRMefu8jTlsg3jluD2zCGTAEXufUyRnfDbOsccAc4Ec+wjo5vwI8DjzSyPLFoIPUPagbWB3oFed7nDo54zXneCCnTs7YSmgnZwQj6zC0mzU32PoQRlk7uY6vAhY5xxnANkfmts5xRqhkDVZeJ14/dNKAhKtsnfvk4H8SwWWcOjljSbjKNQhZu6Fjw2O8/JOBVNfxAnRj3nDKmuV57qgS3emUcVB1J2gZQp3JcDl0fGWO8wf5yFP50K6hZ5zjMcAapxDXALe6ru8JLHEqzOueP10YZb0RqAJWutxQJ+xjR/7PgZeAlBDIOAXYiL7w73X8HkRbLAAJTjltdsqtp+vae53rNgCXNsGzDyTrR8ABVznOClQfwijrb4G1jkyfAP1c137LKe/NwDdDLWsw8jrnD+D18dTUZYu29vY5/5ndaBfbncCdTrgAf3XysQbIC1e5BiHrM8AxV30tdPx7OuW5yqkj9zYDWe9y1ddFuJStr7pzps6afLJYLBZLRBHRY1wWi8ViaX1YxdVKELW4ke8nbLqI/DpE993uTISwWCyWRsEqriAQkRKXqxWRk67zG8ItXzAYYwYaYwrO5FpH+XjyvN9RdCmNLGIwcrQVESMiC738nxSRxxrxPjkiMltEjjn5fcI1w7Qh6RQ4aTT2uromQdQU2UwRKRWRHSJyfT1xS7xcjYj8xUe8XBEpF5GXXH4Fjp/n2g1e1/gNr+++InKXs26sQkSm+5Clv4h8LGqiaLOIXBVseBBp+61DIvKSiOwTkeOiZrFu87q23nKvLzyIPAW6d6PU/ZDTFAO6Lcmhs6Um1RMeE24ZzyBP04FfB5NndNbQKuDhxiivBso5AR0YLgayXP6L0XVtjVUes50ySXDyuwb4QQPTyAFqUCsDXw/RcwtpXUMH4l8FUtC1Q8UEMRPMiV8CjPcR9gHwGfCSy68AuK2e9OoN93df4KuoFZq/A9O9yw6dKPAT1BzRRKAU6BNkuN+0A9UhdKatZ4ZlP3QN5ohgy91feCCZg7z3Wdf9pnC2xdUIOC2Sn4vIaqBURGJE5G4R2SIiJ0QNul7lFf8ex/+YiDwvIgmu8GwR+beocdVtIvIDx/9ar6/LChEpcML6O1+mRaLdglf4kHGSczxMRJY7sr2KVtKgMMbsB94HhrrSNiLS23Xut+vRX96CZChQiJrhudJJLxo4B1jRgHQC0QOdzl/u5Pc99A/fEG5CZ1VNB252BwTx/P2G+6lrPp+9iPQSkaPiGJB1yv6Q+Oky9pIxGbga+IUxpsQYMw+YBXwjiLxfja5F/MwrzWlAETqDNhSccl9jzJvGmLfQ9YTe9ENtfj5mjKkxxnwMzKcuf/WGB0gb6qlDxpi1pm5tnnFcLwhc7gHCA+Wp3nsHkrs5YRVX43EdujYk3ajJoy3A+ehC3F8BL4mzyNjhBtQuWi/UgO59ACIShVpLWIWaRLkQ+JGITDbGvGqMSTHGpKAVdCvwiojEOtd8AHREjZu+LCJ9vYUUtXL/Fmp9IwOdwn51sJkUkS6ooc/NwV7jutZv3oJMYhg6Hfgt6uw59kPr8To/93zHeaH7cu/4uc/jwDQRSRKRzmh+3wtSRg83AS87brKIeJvx8vn8gwz/sq6h07p9PnujVil+jta9JHRx+gvG1WUsIn8Tkb/5kL8PUG2M2ejyW0VwL7GbcewAuu6Thk6b/4mfa34rIodFZL4fxRoo3Od9G4gAg84i3M3j1FOHnHIvA9ajvQiznaBA5d7Q53KazPXcO6DczYZwN/kizeGj68vx+1aA61aiq9498e90hU0BtjjHo4GdXtfeAzzvOo8C3gH+7pyfjzb5o1xxXgEe8JYbte68l1MXhy4gcFdhCXAC/UKbg2OZ3gk3QG/X+XRPepzazRgwbwHKcDWqZNs68qSiL/hljfyM+wPLgGonb9Pd5RXE9ePQdS7tnfP1wI+9ytPn8w+ifpxS14J89rPQLp/VBLke0ZOul9+3gYIA13VHu0h7ePn/Cfi5c/wAp3YVjnaeZTyqfE4AvYINr+++TtivOb2rMBb98PuZc3wxUAm8H0x4fWkHW4fQ7rxx6EdJbDDlXl94sDL7u3dj1P2mcrbF1Xi4DUgiIjeJ7kdTJLon0SDUYoCv+DvQFhToHzDb3TIA/o9TDe8+jP6RPd1s2cAuU2e70JOmLyOW2cAe49RSV9xATDXGpKLbFvTzykuwBJM3n4hOcOgPrDRqKXsJ+jXoaYU1Ck6r8D3gTdQaQXtUUf6uAcncDHxgjDnsnP8Lr+5C/D//YMLdYcE8+3+g9e8vxmXeKgAlQJqXXxqqNOrjG8A8ozsBALqXHPrR5HMCjTFmsTHmhNEdBF5Au7emBBvu7771YYypQlvtl6GK/3/QrYx2BxNeH8HWIaPdefNQ80ffcbwDlbvf8IbI7OvejVT3mwSruBoPd7dId/RlcRfQzhiTjlq1EFd8t92ubmgrCPSltM0Yk+5yqcaYKU7a09Cuoq85FRXn2q5OxXOn6csW2D6gs4iIV9zgMmnMXPQr7A8u7zIgyXXuy5gwBMhbAAY599nqnL+F/kmHUc/4lqiFau+ZZx73ro9LMtDyeMJ5UR5Bu9iCkRERSQSuAS4QnZW1H/gxMERcGwDi//kHE+7+6Kj32YvO/nwceBZ4QEQygskHOsgfIyJuA6lDUKsI9XETddvbeMhHJ6vsdMrjf1GL98v9pGE49b8STLiv+9aLMWa1MeYCY0w7Y8xk6qzlBBVeDw2tQzHUjTMFKvd6w89AZve9z6ruNynhbvJFmsN/V+Ek1/kAoBzoizbJv4k2vW9zxV+Dfu1kAPOA35i6JvxydGwi0TkfBIxEX9KHcEw9ue4Xh77Q70a7CPLRL7R+3jI6cXcCP3TifhXt1gpqVqFz3gGdrTTEOZ8PPOLIeglwEt9dhX7z5oRPx0e3ixN2G/CZ6zwHHeg/ilrIbsxn7CnLGHQcaSbwL1d4fXJe58jUDVXgHvcp8MdAzz+I+uH9LOp99qjC8mzq+TSODckgy2EG2u2YDIwlwKxC1KxTKY79PJd/kldZ/AHds62DU76T0QlCMWjXr3v2Xr3hAe4b41z3W3RMNwHXTExgsOOXhCrTbbi6UusLDyJtn3UIHYechs4IjHbyVsqpJrPqLff6wgPIHMy96637zcWFXYBIcwShuBy/h9EX2GHg/wFzOVVx3QN8gb58X8DZkdcJz3Yq5n7URtkiVOk8gCrAEpd717lmoHOPYifdq/zJiNpAXIG+4F51XNCKy/H7O/BvV3prnfRedGQ/TXHVlzcnbA7wbT8yPIF2dbn9VgK1eL2wGuEZD0XHDI45z+81dM80gpDzPRwF5eV/jZPnmCCev99wP8/C57NHZ17uoc72ZQo6qeYG17VP4hhH9iFzBtqyLUU/dq73Cn8X+D/X+VPAi0GU7wM4Y1yo8lrq1J0ipz5c5Ipbb3h993XuY7zcA67wR51nXOLkpbfX9X7Dg0jbZx1y8jPXyctx9APl2173DVTufsMDyBzMvX3K3Zj/r8Zw1lZhGBCR7agS+yjcsjQXnNmOq4DBpq4LtNnRGHIGev62flgs9dP8VkRbWiVGt/PuH245AhEpclosLRk7OcNisVgsEYXtKrRYLBZLRGFbXBaLxWKJKEI2xiUizwGXAweNMaeZSXHWEf0JXSNQBtxijPG3ruNL2rdvb3Jycs5KttLSUpKTk88qjaYiUmQtKqti//Fy0mNrKKqKJistgfSk2HCLVS+RUrZgZQ0VVtbQ0FiyLlu27LAxpoO3fygnZ0xHpzD/00/4pUCu40aj06tHB0o0JyeHwsLCsxKsoKCA/Pz8s0qjqYgEWd9asYd73lxD+6oa/uecav64JobY2Gju++o5TB3my3hHcCyd9RRdlz9KR3OIg9KBXcN/ysgr7mg0uSOhbD1YWUODlTU0NJasIuLTqk/IFJcx5lMRyaknypXUGcRcJCLpItLJGLMvVDJZ6scYKCmB48fVlZSoKy2FsjKoqIDycnXV1eqqquDJggqOn+wBwLubayk6GEWxGH5YeJKtF0FcXJ1LTFSXlATJyZCaqi4tTV10tMqydNZTDFp2H4lSCQJZHKLNsvtYCo2qvCwWS+QR0skZjuJ6x09X4TvAI0btZSEic1AjnKc1p0TkduB2gMzMzBEzZsw4K7lKSkpISWnyfRDPiLORtbYWiopiOXo0niNH4jh6NI5jx+I4diyW4uJYioriOHEihuPHYzlxIobS0hhqa+uztBNaRAxJSTWkpFSTmnicNiknaZNaTpvUk6SnldG2TRlt2lTQoXcW7dtXkpFRSXT0mdff1lIPmhora2hojbJOmDBhmTEmz9s/ItZxGWOeRs3VkJeXZ862CdpSmtzl5bB9O2zdCjt21Lldu2D3bti7V1tE3iQnQ4cO6nJyICNDXXo6tGmjztMSSk5W52kpxceri4uDmBh1E/7wCXuKTwKGn5xTzf9bEwu1QkZ8Ir8cmU9lpcpRVaWtNo8rL4eTJ9WVlUFpqVBSEkNJSQwJOz/hSFE7du5tx+GyHIrK00/LhwhkZkK3btC1q/527655ysmBHj20FefNfwt+wZ+2zmRqm+/y1sG/8cOeV3FZ/kNn+oiahJZSZ5sbVtbQEGpZw6m49nCqIdEu+DYK26qpqhK++AI2bICNG2HTJnVbtsAer9KKja17iY8fD507q8vOhk6dICtLX/RJSb7vdab8bEof7nlzDSeraoiKAokyxMdGcdOIXLp2DXy9L7Lf/j6dOPTleUV1HAdLO7CmuB9Lur3J0aNw5AgcPgxHj8KyZTB7tipBN+3bQ+/e0KsX9OkDJdWv8U7UeuikhbAvWnhg20yAZq+8LBaLEk7FNQu4S0RmoJMyilvz+FZFBaxbB2vX1rl162DLlvHUujas6NgRcnNh0iR9GffsqS2LnBxVTFFhWODgmYDxyOwNwAnaJyVyw6C+jO9+5hMzVnT9Kek7nTEuID6mknaphzF9pvCVUafG9Yy1VVSoMtu3Dw4e1N99++DAAZgzB/71LzDmGtRsINzX/gRVWYOJz97Gz7I30T4RBgzQlqbFYmm+hHI6/Cuoper2IrIb+CVqvRpjzJPorptTUKOfZagF9RaPMbB/P6xcqW7VKlizRltUNTUaJyZGldM558Do0TuZPLk7fftqi6FNm3BK75+pwzoztktnFi4s4KnL8s86vcyhd/AJMGzXo2Saw+yv6cCKzj8na9Rtp8X1dFkmJmp3Z69ep4bX1KhiKymB7y+5ivIDOVTs60H3Q5ewZmcGpQUjOFKZyLnTNX5ODgwZAoMH6+/QofqBIOEb/rNYLC5COavwugDhBvheqO7fHDBGx5sKC2H58jp34EBdnJwcfUFedZUqqkGDVGnFxWl4QcE28vO7h0X+cJM59A72Dr2DV+aq8sg6Q6UdHa0uIQG6Z29if7fNANySEsPfS/6OqRUS9/ZgSuV/2L5dxw1XroS33+bL1m5amsowfDjk5anLza1r4RY8+yCxT79GenENRW2iqbr9GvJvvf/sCsBisfgkIiZnRApHj8KSJXVu6VLtsgJ9cQ4YAJdeCsOG6Utw8GBtIViajutTr+JvpTMpj6prPsUbw03dz2OMy3RuZSUUF9eNJ27erJNg/v537ZIE7VLMy4Os+PmMWHuYEbHtiIo9SEZxDRWPv0IBWOVlsYQAq7jOkJoaHYdauBAWLIBFi3TyBGiXUv/+qqTy8mDkSFVSiYnhldkCg3s9xHe3wEvH/wNAu8oYrk68ljH97z4lXlxc3czLMWPUr7pauxs3bdJn7Zks89mmUbxixgKQFbOPwQmrGJK4il5/XM15N+osTIvF0nhYxRUkpaWweDHMmwfz56vCOnFCwzp2hHPPhVtugdGjVVn5moZtaR4M7vUQ/3PgIUpLC/jj8BVBXxcToy3kkSPVgXYlRt+Qx8aKfqwuH8zqk0NYVT6ED0ougUPwwzYwYgSMHQvnn69KsF270OTLYmktWMXlh2PH4LPP1H36qY5NVVdra+qcc+DGG/UldN55duC+NRMVBWVtaxhSvJohiauh7UsAHKpuz6e1o1g5/I988QU89hg8+qhe07+/LlcYPx4uuECXLFgsluCxisvhyBGYOxcKCvR3zRqdXBEXB6NGwU9/ql/M551nx6Usp7J36jWkvPwK8dV1fqnxR8n9WicuvULrUXGxTvj4/HN1L70ETz2lcXv0gAkTVInl5+taPA8rn5zNsqXlpF0Qy7O3vsmIkQkMvXNKU2bPYml2tFrFVVysLamPP1a3erX6JyZqS+pXv9IXyahROhvNYvFHr4vvZyOQ/dZrtC2u4UhSAvu+cgO9r/hfQFvj6emqlPLzVZGVlakC8yyHeP11eO45TS8nBy68EHrWrkTKE0lLSyCNUspj01m4rBKenG2Vl6VV0+oU13/+A3ffPZyNG3V8IiFBFdVDD+lX78iRdVPRLZZg6XXx/XDx/Xy6RFtMvbP8xxVRM1qjR6sDtfjxxRfaJb1qFbz6KpSUDAUgK72UIZsP0i36BLnZRSxbepihd4Y+TxZLc6XVKS5jICbGcO+9+lU7erRtUVnCT2KiTuIYMULPy8th1h+XsHFfWzbuyaDgs65UVMYgYujW/jhb71XrKWPG2FmLltZHq1NcU6dCevqKiDFWaWmdJCRAbvpuunUsYdKQXWSMKmPJfzqxYU8GG3a14Xe/g9/8Ru1Onn8+TJ4MF1+sawXtRCFLS6fVKS6LJVLolpPAll2V1EbHERNj6N2pmN4dD3PH5XEkjpnAokVqXHj5cnj/fb2mU6c6JXbRRWpk2GJpaVjFZbE0U9pdNgX+O5sd24uBGKIrSunUPYUuX5kAwGWXqauthW3bdG1hYSH8+98wfbq2vIYPh0suUXfuuboW7ePpb7Hqw9eorT5OVEwaQy66hom3TA1nVi2WBmEVl8XSjGl32RTiS2D//gJG3pXvM05UlBoW7tVL1xdWVGgrbPFibZE98gg8/LAuih+Su4cs2U3frFjSk6C2+jgr3p0OYJWXJWKwistiaWHEx+t6w/PO0/NDh9Tay9KlsOizBD4ruwuArDbb6Zu1jH5Zy6l9702ruCwRg1VcFksLp0MHnZQ0dSp88rubOXC8G+v3j2DD/uHM23QFczdcTVzMST67AqZMURub3VvnhgSWCMEqLoulFREVk0ZWm51ktdlJft+ZVFQlsPnQYNbtG8OSJZN4+22N17evjp9NmQLjxtkp95bmRRj2y7VYLOEie9g1uL9X42PLGZC9gju/Vcbbb8M//wm3367T7P/yF10r1q4dXHmlmqjauTN8slssHmyLy2JpReROnArA3hWvUVtznBqTTtbQaxg0+QpA14ENGKDKq7hYx8YWL9Zte2bN0jT6969rjY0dW2dpZtOMAszyk1QNqWHjz95FhieSOy2/6TNpafFYxWWxtDJyJ04ld+JU1qzRllWvXr7jtWmjymnKFJ1yv26d7pawdCk8/jj84Q+QkqLrxUa038CF0TF0aZMCFJMUlUL18io2UWCVl6XRsYrLYrEEJCoKBg5Ud+edUFSke9MtXqzKbObhvtxHX/p1KGXI6H1cnVpGXufjmOUnYFq4pbe0NKzislgsDSY9HS6/XF1tLRT/vZCCbRl8sjWDN97tyas1UaTGVTMu5xhf7aMLoLt0CbfUlpZCSBWXiFwC/AmIBp4xxjziFX4L8Ciwx/F6whjzTChlslgsjUtUFHRrd4g7O57kztF7KOxVyqF3u1GwtS1ztrTl3W9rvAED6roe3WNjFktDCZniEpFo4K/ARcBuYKmIzDLGfOEV9VVjzF2hksNisYSew90Syd5VRUxULMlJ1eT1OcKk3gf4VmYCu9qdx7x5ao7KMzaWnKzbCE2ZorYVe/YMdw4skURAxSUi3wdeMsYca2Dao4DNxpitTjozgCsBb8VlsVginPgx+exdUED7nSUAnKg6yf7OqaRfcB7pwDnnaDz3TMWlS+Gdd9S/Vy/tTrz0Ut1sMzm5Lu2ls+Yxd/k8Skw5KZLABcPHMfKKcU2ZPUszQ4wx9UcQ+TU6vLoceA543wS6SK/7GnCJMeY25/wbwGh368rpKvwtcAjYCPzYGLPLR1q3A7cDZGZmjpgxY0ZQmfNHSUkJKSkpZ5VGUxFJslZXQ2lpCbGxjSvviRM6+y06uvHSrK6GmpoS4uMbV9bSUu0Ci41tvDRra6GqqvFlPXlSu/kae3FxRUVwshoDu3cnUliYwYoVGaxZk05lZTSxsbUMGlTMyJFHGdRvDxkd9hDlfvYG2iSnktTm7Msjkv5frVHWCRMmLDPG5Hn7B1RcACIiwMXAN4E84DXgWWPMlnquCUZxtQNKjDEVInIHcK0xZmJ9suTl5ZnCwsKAMtdHQUFBxOzHFUmyHjoECxcW0LlzfqOmO3cuDB2q07MbiwMHoLS0gJ498xsvUWCJswNyVj07IDeUEsfIbu/e+Y2XKAScDn+mbNpUQG5ufoOvKy3V8vNs17J9u/onJ5fQs+dWevbcQq9eW0lNPUEKCfzvA3eftayR9P9qjbKKiE/FFdQYlzHGiMh+YD9QDbQF3hCRD40xP/Nz2R6gq+u8C3WTMDzpHnGdPgP8Phh5LBZLy8Mz7jVBd21h926Y/vxMtmztxZYtvVizZjAAHTocpGfPLQwYBePH61oyS+simDGuHwI3AYdR5fJTY0yViEQBmwB/imspkCsiPVCFNQ243ivtTsaYfc7pFcC6M8qFxWJpcXTpAmOGbmDI0FUYIxw4kMmWLT3ZurUXywpHctllur/Y6NG6ceaFF8KoUY3bTWtpngTT4soAvmqM2eH2NMbUisjl/i4yxlSLyF3A++h0+OeMMWtF5EGg0BgzC/iBiFyBtuKOArecYT4sFksLZGDXcSzb+TE1UktW1n6ysvZz/nmL6dH2Yg7UjGbJElixAh54AH75S221jR+vNhYnToTBg3Ucz9KyCKi4jDG/rCes3haSMWY2MNvL737X8T3APYHFtFgsrZHMoeMYAazeMZ9yOUlcbRK52efTc9RoBlDXrXj4MCxYoFPuV62Cd99V/4yMOkU2YYLaWRSB1atXM2fOHLKysnjssce48MILGTx4cLiyaWkg1nKGxWJp1mQOHUfP1HGUldVNq/emfXu44gp1ALt2wcKFuhP0ggXw1lvq36EDDBtWRHz8Hrp2jSMzE4qLi3nb2c/FKq/IwCoui8XS4ujaVd011+i0+61bVZGtWgWLFkVz/PilALz0UjmdO2+le/ftHDu2it//frDtWowArOKyWCwtGhGd8t+rF9x4I8ya9UeOHWvL9u05HDs2mDVrslm3bgDvvQfPPqvmqM4/X92IEXYTzeaIVVwWi6VVkZTUBpFjZGQco2/fMjZs2EBRUTqHDvUjOfkSPvsM/vtfjRsfrzMVx46FMWPUtWsXXvktVnFZLJZWRr9+F7J69dvU1FR96dehQynf+lY2niGuQ4fUNNW8ebptyx/+oNZWAPr2hfPOUyV23nlqPNh2LzYtVnFZLJZWRZcuqp3Wr58DQFpaGyZNOnVWYYcOMHWqOoCyMp2xOH++TvZ4+22YPh3nem2VjR4N556rxx07Nl1+WiNWcVksllZHly6D6dJlMHv2FPC9710XcBwrKUmn1Y8fr+fGwObNOuFj0SI1GvzII1BTo+E5OarARo5UN3w4pKaGNEutCqu4LBaLpYGIQG6uuptuUr+yMrWxuGSJKrLFi+G11+ri9+0LeXk64WPECLW/aZXZmWEVl8VisTQCSUl1sxE9HDyoXYyFhbqNy5w58NJLGiYCffrAsGF1buhQ7aa01I9VXBaLxRIiOnas2/XZw7592jJbvlzNVS1YAO6dmrKzYcgQdYMHq+vTp+llb85YxWWxWCxNSKdOcPnl6jwcPQorV9a5Vavgo4+gypn4GBcHXbvmcd55MGiQuoEDdQud1jij0Soui8ViCTMZGWoUeKJrN8LKSli3TvdNW70a5s6t4JNPUr7sagQ1Kty/v07JHzBAj/v1g5491XJ+S6UFZ81isVgil7i4ui5DgIKCNeTn51NUBGvXwuefwxdfqPvoI/jnP+uujY3ViSN9+9a5Pn3UtWun42uRjFVcFovFEkGkp6slj7FjT/UvKoING7SVtm4drF+vv2+/Xbd42nN9bi707l3nevXSVlpWVmQoNau4LBaLpQWQnq6LoEePPtW/uhq2bYONG2HTJv3dvFnXn736KtTW1sVNSoIePU51OTl1Lj29eSg2q7gsFoulBRMTU7fmzJvKSlVqW7eq27Kl7nzuXDhx4tT4qanQvbtOCunWrc4Kv8d17gwJCU2Qp9DfwmKxWCzNkbi4ujEwb4zR2Y7bt8OOHarQdu5Ut2OHLrA+cuT069q1g/T0PN5917eybAys4rJYLBbLaYioEmrXTi19+KKsDHbv1o07d++uc6tWldOmTUrIZLOKy2KxWCxnRFJS3WxFNwUFn9OxY37I7tsKl65ZLBaLJZIJqeISkUtEZIOIbBaRu32Ex4vIq074YhHJCaU8FovFYol8Qqa4RCQa+CtwKTAAuE5EBnhFuxU4ZozpDTwG/C5U8lgsFoulZRDKFtcoYLMxZqsxphKYAVzpFedK4AXn+A3gQpHmsErAYrFYLM0VMcaEJmGRrwGXGGNuc86/AYw2xtzlivO5E2e3c77FiXPYK63bgdsBMjMzR8xwm1I+A0pKSkhJCd2Ml8YkkmStqVF5RRpX3spKNWHTmJ80tbVgTAnR0Y0ra1UVREc3ruFTY6C2tvFlra7WMo2ObtRkqalpfFlrarQcGtv+njElpKWlNItFtYGIpHdBY8k6YcKEZcaYPG//iJhVaIx5GngaIC8vz+Tn559VegUFBZxtGk1FJMkKKu8FF+Q3aprGhGa1/ty5jS9rqLCyhoa5cyPn/xVJ74JQyxpKxbUH6Oo67+L4+YqzW0RigDaAjyVtlkiisZVMKL+GI+FL24OV1WJRQqm4lgK5ItIDVVDTgOu94swCbgYWAl8DPjYB+i6XLVt2WER2nKVs7YHDAWM1DyJJVogsea2socHKGhpao6zdfXmGTHEZY6pF5C7gfSAaeM4Ys1ZEHgQKjTGzgGeBF0VkM3AUVW6B0j3rja1FpNBXv2lzJJJkhciS18oaGqysocHKWkdIx7iMMbOB2V5+97uOy4Gvh1IGi8VisbQsrOUMi8VisUQUrVVxPR1uARpAJMkKkSWvlTU0WFlDg5XVIWTruCwWi8ViCQWttcVlsVgslgjFKi6LxWKxRBQtVnGJSIaIfCgim5zftn7i1YjISsfNcvn3cCzWb3Ys2MeFU1YRGSoiC0VkrYisFpFrXWHTRWSbKx9DQyDjGVv6F5F7HP8NIjK5sWU7A1l/IiJfOOU4R0S6u8J81ocwynqLiBxyyXSbK+xmp85sEpGbQy1rkPI+5pJ1o4gUucKarGxF5DkROeiYlfMVLiLyZycfq0VkuCusScs1CFlvcGRcIyILRGSIK2y7479SRAqbgaz5IlLses73u8LqrTsNwhjTIh3we+Bu5/hu4Hd+4pX48X8NmOYcPwl8J5yyAn2AXOc4G9gHpDvn04GvhVC+aGAL0BOIA1YBA7zifBd40jmeBrzqHA9w4scDPZx0osMs6wQgyTn+jkfW+upDGGW9BXjCx7UZwFbnt61z3Dbc8nrF/z66fjMcZTseGA587id8CvAuIMC5wOIwlmsgWcd4ZEB321jsCtsOtG9G5ZoPvHO2dSeQa7EtLk61PP8CMDXYC0VEgImoxfoGX38GBJTVGLPRGLPJOd4LHATOejF2kJyNpf8rgRnGmApjzDZgs5Ne2GQ1xnxijClzTheh5sjCQTDl6o/JwIfGmKPGmGPAh8AlIZLTQ0PlvQ54JcQy+cQY8ylq1MAfVwL/NMoiIF1EOhGGcg0kqzFmgSMLhLe+BlOu/jibun4aLVlxZRpj9jnH+4FMP/ESRKRQRBaJyFTHrx1QZIypds53A51DJ2rQsgIgIqPQr5YtLu+Hne6Ex0QkvpHl6wzscp37Ko8v4zjlVoyWYzDXNiYNvd+t6Je3B1/1IVQEK+vVzrN9Q0Q89j+bulwbdE+n+7UH8LHLuynLNhD+8hKOcm0I3vXVAB+IyDLRXTSaA+eJyCoReVdEBjp+jVquEWEd3h8i8hGQ5SPoXveJMcaIiL95/92NMXtEpCfwsYisQV+6jUojyYrzVfgicLMxptbxvgdVeHHo+omfAw82htwtGRG5EcgDLnB5n1YfjDFbfKfQJLwNvGKMqRCRO9BW7cQwyhMs04A3jDE1Lr/mVrYRhYhMQBXXOJf3OKdMOwIfish6p1UULpajz7lERKYAbwG5jX2TiG5xGWMmGWMG+XD/AQ44L3nPy/6gnzT2OL9bgQJgGGqhPl3UYj34tmzf5LKKSBrwX+Bep3vDk/Y+p8ujAniexu+Ka4ilf+RUS//BXNuYBHU/EZmEfjRc4ZQb4Lc+hE1WY8wRl3zPACOCvTYENOSe0/DqJmzisg2Ev7yEo1wDIiKD0ed/pTHmyx00XGV6EJhJaLvhA2KMOW6MKXGOZwOxItKexi7XMx0ca+4OeJRTJzz83kectkC8c9we2IQzYAi8zqmTM74bZlnjgDnAj3yEdXJ+BXgceKSR5YtBB6l7UDewOtArzvc4dXLGa87xQE6dnLGV0E7OCEbWYWg3a26w9SGMsnZyHV8FLHKOM4BtjsxtneOMUMkarLxOvH7opAEJV9k698nB/ySCyzh1csaScJVrELJ2Q8eGx3j5JwOpruMF6Ma84ZQ1y/PcUSW60ynjoOpO0DKEOpPhcuj4yhznD/KRp/KhXUPPOMdjgDVOIa4BbnVd3xNY4lSY1z1/ujDKeiNQBax0uaFO2MeO/J8DLwEpIZBxCrARfeHf6/g9iLZYABKcctrslFtP17X3OtdtAC5tgmcfSNaPgAOucpwVqD6EUdbfAmsdmT4B+rmu/ZZT3puBb4Za1mDkdc4fwOvjqanLFm3t7XP+M7vRLrY7gTudcAH+6uRjDZAXrnINQtZngGOu+lro+Pd0ynOVU0fubQay3uWqr4twKVtfdedMnTX5ZLFYLJaIIqLHuCwWi8XS+rCKq4UiamEj30/YdBH5dYjuu92Z+GCxWCwhwSouH4hIicvVishJ1/kN4ZYvGIwxA40xBWdyraN8PHne7yi6lEYWMRg52oqIEZGFXv5PishjjXifu5z1RRUiMt1PnGkisk5ESkVki4ic38B7FIjIsRCssWsyRE2TzXTKYIeIXF9P3BIvVyMif/GKkysi5SLyko/r/Za3U5blrrQ3BBPmhPt91kFcmyMis53nuF9EnnDNPEZE+ovIx6ImjzaLyFXB3DfIPAcse3/lGejeAcrEb57CiVVcPjDGpHgcOivmKy6/l91x3RW3hfEVJ/9D0Vl494RBhqHo+rQBIuJeAzcMHaRuLPYCvwae8xUoIhcBvwO+CaSiZm+2Bpu4qN3G89HFolecpaz13SfUdfGvQCW6QP4G4O9St8D0FLz+Q1nASXTyjnd6S72vDbK873Ldo28Dwup91gGu/Ru6VKUTWjcvQE2decr+P8A76MzE24GXRKRPMPcNIs/BlL3P8gwizz7Dg8hT2LCK6wxwWiQ/F5HVQKmIxIjI3c5X0glRA65XecW/x/E/JiLPi0iCKzxbRP4takx1m4j8wPG/1uurtUJECpyw/s4XYpFot+AVPmSc5BwPE5HljmyvojMAg8IYsx94H/2jetI2ItLbde6369Ff3oJkKFCImt250kkvGjgHWNGAdOrFGPOmMeYtdN2ZL34FPGiMWWSMqTXG7DHO+pkguQmdYTUduNk7sL76EUTd8VUXfdYNEeklIkfFMSjrPJtD4qdL2UvGZOBq4BfGmBJjzDxgFvCNIPJ/NfrC/8yV3jSgCJ1N683ZlrdfgnjW9dEDXeZR7vwv3kOXe4AuAcgGHjPG1BhjPgbm45TP2dSxYMq+vvIMdO96wuvNUzixiuvMuQ5dC5Ju1MTRFvSrug1aCV8SZ1Gxww2oHbReqMHc+wBEJAq1jrAKNYFyIfAjEZlsjHnV9dWajX6BvSIisc41HwAdUWOmL4uI9xciolbt30KtbWSgX71XB5tJEemCGvbcHOw1rmv95i3IJDwtq7eos9/YD6236/zc8x3nhe3LvXMGeYhGlyV0cLpKdot2ESU2IJmbgJcdN1lEfJn08lk/gggDV11Ep3n7rBtGrVT8HK2bSehi9Rc8Xcoi8jcR+ZufPPQBqo0xG11+q6h7cdfHzTh2AZ37pKFT6H/iHbEB5f1bETksIvN9KN76wgJR37WPA9NEJElEOqP/i/fqSUuAQYFuGESe6y37+sozBASVp5AT6nn/ke7QhZSTfPh9K8B1K9FV7p74d7rCpgBbnOPRwE6va+8BnnedR6HN9b875+ejXWhRrjivAA94y412Oezl1MWgC4BfB8hzCXAC7d6ag2OJ3gk3QG/X+XRPeu7yCiZvAcpwNapk2zrypKIv8WUheta/BqZ7+WU7+S1Eu4jao1+dDweZ5jh0zUt753w98GMf5e2vfvgN81UXg6wbs9C1S6sJcn2iJ10vv28DBQGu6w7UAD1cfn8Cfu4cPwC81JDydupVKrqo/WannvYKFBbEs673WqA/sAyodmScTt1i21j0w/JnzvHFaNfe+2dbxwKVfX3lGeje9YUHm6dwONviOnPcBiMRkZtE958pEt2DaBBaAX3F34FWVtA/dra7ZQD8H6ca2n0Y/UN5utmygV2mzlahJ01fRiuzgT3GqYmuuIGYaoxJRbcp6OeVl2AJJm8+EZ3E0B9YadQy9hL0C7exx7cCcdL5/YtR01qHgf+HKpBguBn4wLkO4F/46C7Ef/0IFOYdHkzd+AdaP/9iXOauAlACpHn5paEv9/r4BjDP6M4AiO4VNwnwN7kmYHkbYxYbY04Y3XHgBfQlPyVQWCDqu9bpPXgPeBO1UtEe/aD6nXNtFdorcBn64fA/6NZIu4O4daA8+y37IMrzjDnLPIWUljqxoCn4UhGIWsL+B9oVttAYUyMiK9FmtQe3na5uaCsI9KWzzRjj0xCl03d9HTDSqUg413YVkSjXC6obuirdm31AZxERl/LqxqmW5f1n0pi5ojON/kBdd10ZkOSKloXvylxv3gIwyLmPZ4D6Lef+mahNNp+IyLvoF6ovPjPGXNoQIYwxx0RkN67n7XXsF6er5xogWkT2O97xqB3MIcaYVa7o/upHoDBveeqtG6KzQx8HngUeEJF/G2OC2aZiIxAjIrnG2V4HGIJaSaiPm4BHXOf5qMmgnSICkIKWzwBjzPAzLG/Dqf+1YMMC4b42Ay3HJxxlXyEiz6OtlJ8BGGNW4zLYLCILqNvqx/9NAue5vrLPp57ybFBufct2RnkKOeFu8jV3h/+uwkmu8wFAOdAX3TDtm2h3wm2u+GtQw5IZwDzgN05YNGpR+edAonM+CBiJti4O4Zh2ct0vDn2h34024fPRL99+3jI6cXcCP3TifhXtugrUVejOXwegFBjinM9HX0bR6F5FJ/HdVeg3b074dPx0XQC3oYrGc56DDj4fRS1iN+YzjkEnrPwWHQtMAGJc4Q+is7U6ol/ZnwEPucJ95gP94DiKvvCyXO5T4I9e5e2vfvgN8/Os6q0bqMLybPL5NI5NySDLaQba7ZgMjEV3UfBrbw4181SKY0/P8UvyKos/oPu3dQimvNFxvMmeZ4R2HZei40B+wwI96yCv9ZSrJ/5M4F+u8MHO9UnA/6J2DuPru28D6pjPsg+yPAPd2294fXkKpwvrzSPBEYTicvweRl9Snmb+XE5VXPcAX6Av3xdwduB1wrOdSrkftUm2CFU6D6AKsMTl3nWuGejco9hJ9yp/MqIDvyvQF9irjgtacTl+fwf+7UpvrZPei47spymu+vLmhM0Bvu1HhifQrhO330qgFteLsJGe8QPoF67bPeAKj0WnQhc5+fgzkOAK95kPtGvpjz78r3HS8bwc/NaPIOqOr2fls26gMzP3UGcLMwWddHODc/4kjqFkP+WUgbZ8S9GPoeu9wt8F/s91/hTwYhBl/5KXn9/yRj+iljp1r8ipTxcFCgv0rIO8dihq1f4Y+j9/Dd1LzxP+qBNW4pRF70D3bUAdq7fsA5RnoHv7Da8vT+F01lZhEyAi21El9lG4ZWkuOLMdVwGDTV0XaMTRGPmor37YumOxnI4d47KEBaPbd/cPtxxnS0vJh8USSdhZhRaLxWKJKGxXocVisVgiCtvislgsFktEEbIxLhF5DrgcOGiMOc1EiOiigz+hi+zKgFuMMcsDpdu+fXuTk5NzVrKVlpaSnJx8Vmk0FZEia1FZFfuPl5MeW0NRVTRZaQmkJ8WeFq+06BBxZQeIpYoqYqlMyiQ5vUMYJI6csgUra6iwsoaGxpJ12bJlh40xp70gQjk5Yzo6pfmffsIvBXIdNxqdbj06UKI5OTkUFhaelWAFBQXk5+efVRpNRSTI+taKPdzz5hraV9XwP+dU88c1McTGRnPfV89h6rA6gw1LZz3FoGX3kSjx6DpcOGnK+HzEtxl5xR2npfvfgl/wp60z2R8FWbXww55XcVn+Q40mdySUrQcra2iwsoaGxpJVRHxa+QmZ4jLGfCq6nYM/rqTO8OYiEUkXkU7GmH2hksniG2OgrAyKiqC4GI4fV1dSAqWl+nvyJFRUQHm5uupqqKpSN3M5lFYMAAOvLK7lyDGBKPjOnFo+PRfi4tRVzCvl05jbSYwtJym2jOTYMlLjT1C5aw6rc+4gPR3atoWUFJg99xc8sG0m5dFquGBfNDywbSaAT+VV8OyDxD79GunFNRS1iabq9mvIv/X+JixFi8XSVIR0coajuN7x01X4DvCIURP9iMgc1FDkac0pEbkd3QuGzMzMETNmzDgruUpKSkhJafJ9Ec+IM5W1pgaKiuI4fDieI0fiOHIkjqNH4zh2LI6ioliKiuIoLo7l+PEYTpyIpaoq+OHO2NhaYmNriY42REcbaj2WhQSiBWpqodYItbUCtVFUVUU1KP2oKENSSjnJqeWkpJ4kJa2clDYnSWtTRnqbk/TrlkH79pW0a1dBRkYlFcV7iT54FHFVZSNQ0zGDlPad/N6nNdSDcGBlDQ2tUdYJEyYsM8bkeftHxDouY8zTqHka8vLyzNk2QVtCk7uqCnbuhK1b1e3YUed27YK9e7VV5E3bttCxI3ToAD17Qvv2kJGhLj0d2rRRl5qqLZ/kZHWJieri4kDtjdYporGPfMyeIrUT6ukqBGiXmMjjEydiDNTWQo85/WhbXcLJqgTKqpIorUriREUK20525Yuez1NaCidOQEmJ8OHB/3CypC0lJelU721L9fosak60BXOqAhSBjJiDZEfvo1PMPrJj99I5dg+dY/eS0nYvFy75D4leG2KsfHI2y5aWk3ZBLAdmHmXEyASG3hmszdzw0BLqbHPEyhoaQi1rOBXXHk41HtrF8bM41NbCgQPxfPABrF8PGzfCpk3qdu7UVpWHmBjo0gW6d4fx46FrVz3v3BmysyErCzIzVfE0Nj+d3Jd73lzDyao6geKjo7nxnL4kuUzxrs35MRN23kdafJ1B8bLaeI52msZVo05N88jKX7M/5lTbqKY2ipQjHbkl5UMOHoSjR+HwYSh75zMOVGWxqbIPc0vzqTDOPot74KYkzX/v3tCnD6QcWUd5SSYdMqoYUHOY8th0Fi6rhCdnN3vlZbFYlHAqrlnAXSIyA52UUdxax7eMUUW0ejWsXQtffKG/69dDWdl5X8ZLTYXcXBg1Cq6/Hnr10lZTjx6qoKKjwyO/ZwLGI7M3ACdon5TIDYP6Mr77qbusZA69g0+AYbseJdMcZl9NR1Z2/hlZo247Lc3rU6/ib6UzKY+qU17xppYbsi5kpJedipObfkm7YlWaxsDhmvbsrerMWunJjvG/Ye9e2LdPy7eoqO7i6Ddq6ZBWRlZ6Ge8WHuGm9jBokCq5GOef8fH0t1j14WvUVh8nKiaNIRddw8Rbpp51mVksljMnlNPhX0EtU7d3TPb/EjUkiTHmSWA2OhV+Mzod/puhkqU5UVYGn38Oq1bBypX6u2aNTobw0KULDBgA3/42REVt4Mor+9K3r7aY5Ew3aAgxU4d1ZmyXzixcWMBTl+X7jZc59A72Dr2D/67Q7snu3X3HG9zrIb67Bf51Yib7o4WMyli+lnQNY/rffVrcvVOvIeXlV4iv1vLpEHOY1PijyNfymHqFxjFGu04/e3IuB48ncaAomRNJsWz5PJ09R1JYtb0Db39d48bHa/l3TN5BfMkRstO70zl9K4kcZ8W70wGs8rJYwkgoZxVeFyDcAN8L1f2bAydPwooVUFgIy5bB8uXamqp15jKkpcGQIfCNb8A556gbOFDHmDwUFOzjggv6hicDYWZwr4cYzEPMnQtDh55aLm56XXw/G4Hsma/R9ngNR5Pj2Xv5jfS+4n+/jCMCsbHQNuoISR1ryOl4go5jSjnYRdeaVJeWEz92Clu2wLZtOm44f2k6JRW3fplGu+R9dG67hU827qK2C4wYoWOGAJtmFGCWnyRRkjlpSpHhieROyw9NwVgsrZyImJwRCdTUqFJavBiWLFH3+ed141CZmfqiu+oqGDZMX8Q5Oc23BRVp9Lr4firz72dWAVx8MfT2E69bTgJbdlVSG+0a7KupIqdXIl1GwsiRdd4Fv7+B4yfT2VvUkz3HerGnqBd7jvVi9e5x/Pci5769oH/WAYbHZzMiu4yBmSUkxaRQvbyKTRRY5WWxhACruM6QI0dg0SJYuBAWLIClS3W9E+hX+KhR8JWvQF6evgyzvTdbt4SFdpdNgf/OZse240A00RWldOqeQpevTDgtrkSnkZZYRFricvp1qjPqcqKiMwnnPsWGDTphZumqNrxTkglAXHQtAzuWMDz7BIPWH+Rr47Tr183SWfOYu3weJaacFEngguHjGHnFuFBm22JpUVjFFQTGaNfRvHnq5s+Hdes0LDpau/tuvhnOPRdGj9bBfduSar60u2wK1QegtLSAkXfl+42XPewa9hROR/fyVAyx9Bh1GYMugoucVlfH15awvySelftSWbE3leV703h5VRblyzrz47d1hufYsTBuHLStWc7mQ59QG10DAiWU8/6yjwGs8rJYgsQqLh/U1mq339y58Nln8OmnOisNdK3TmDFw4436O3KkrnOytDxyJ04FYPey16G2GBPVho6Dr2XQ5CtOiXeytpROqUKn1CNc2ucIAFU1wtI9sXzSZjRr18KHH4Kumx9OfPwAunXbSbduO8jJ2UGnTnuZu3yeVVwWS5BYxYUqqs8/h4ICdZ9+ql2BoN08+flw/vnqBgyAKGtTv9WQO3EqtZ2nYgz06+c7zuFuiWTvqiImymVUWKrokBfNbRfoaXW1Tvp49dU32bmzOzt3dmPTpj4AxMZW0rXrLirjtK7l5dWtt1uy5C8cOfos1VW38e57P6Jdxq2MGvX90GXYYokAWqXiMkbXSH38sbqCgjpF1aOHjk1dcIE6O4HCEoj4MfnsXVBA+50lJEYlc6KqggNdUki7oG4NXkyMrsEbM2QjQ4asBqCkJJmdO7uzY0d3dmzvwb33atzkZO1azM2dT99+BfTtW0p1FcTFnaCo+C8sWYJVXpZWTatTXC+9BD/60XlfKqpu3eDyy2HCBP3a9beuyGKpj/gx+ZwYA3OWaJ3KyvIdb2DXcSzb+TE1UktKSikDBnzBoH4b6Jc1iaRuHVm8WNf2rV4NH3wwFhhLcvJxBg06TF7eUYYPn0919bOAVVyW1kurU1xZWTBkSBHXXZfJxInawrItKktTkTl0HCOAtbvmUWrKiatNIjf7fHqO0tbZ1KnqjIFt20exauV5LF8xhhXLL2HxYrV2n9HuAJdN0ckhkyZBJ8eO8LMbPuSxvdEcNum0lyJ+nF3DrX0vCks+LZZQ0uoU16RJEBOzjvz8zHCLYmmlZA4dR+bQcXz0kc40TEg4PY4IJCdXMmHiO0yY+A5lpYcoPv46K5aPYcnSfGbNmsyLL2rcgQMhO28HK/p3JWpwBZJgOEwGv9pbAXxolZelxdHqFJfFEimY2lupqfkL0dG6ir1Tp910mPwm48bl0KmTmgpbtEgtsnz0ry6YqmiIrSXunFLi8o5TNfIE/8+UcmvrNLxiacFYxWWxNFOysr7P/v1QJc8CUFbWhuqqO+jWTXeLHj5cHcCUA6up/DyNysJUKpamUvJ0Z3gaitr14ubL4NJLtWuxXTv4v4INzJi7jZqTNUQnRjPtgh78Jt9qN0vkYBWXxdKMycr6PiUl32f//gJ6917uN16HxGMcHinEjzxB6neg5lAslYWp1C5OZObMDvzzn9r9mNX3JGXtokjskUpcpyJqT9bwrw+2AFjlZYkYrOKyWFoAN8bU8LeqCiolHoDoDlWkTd7Pty6q4pKUDqxYoabJ3phbQcWG3hTPzyUqoZKEnMMk9jzES5V7rOKyRAxWcVksLYCLki6Csg/5Z2UMRyWd9OoTXCvRXJExFlBTZKNHw2dd5lN7Mpby7e05ubUD5ds6ULY+myOzYehHcJnTrXjuubr2bOmsp+i6/FE6mkMclA7sGv5TRl5xR5hza2ntWMVlsbQQLkq6iH5H4NChurEvb6ISohGqSO6/j+T++zAGqg6mcXJbFtWHc/nd7+A3v9EtZEYO2sI17dfQrTdEpUAWh2iz7D6WglVelrBiFZfF0oqYNKAHH6zcgtQawNmnLOsEU/Kz+M4gOHxYDUkvXgyrFqbyUckT3A6M6LSCS3t/yJTcD+hm/ghWcVnCiLW6Z7G0Iu7s1ZeLh/ZCEqIxAFGxXNQvl+8MygWgfXtdAP3b38KeH/dh+e3n8/DEB0mIKec38/6HMc99xOA/LGLaNHjxRTh4sC7t/xb8goufG8QXR77g4ucG8d+CX4Qji5ZWgG1xWSytjDt79eXOXn1ZswaSknQzTF8cimrPsE6rGdZpNf93/h85djKdD7dO4I2NV/D++1/l1Ve1xTZ8OPQZ+Akrum4mulc0APuihQe2zQTgsvyHmiprllaCVVwWi8UnK7r+lPSd95EolQC0TSzi8v6zSZ44jv8ZrjYV582DwkJ45cXxYCYQnXKM54ceomjALlIGzedPW2daxWVpdKzislgsPskcegefAMN2PUqmOcy+mg6s7PxzskbdBsCIEeoAvrFwPCfWjqVk9TjWr55EybzfgtSyI2ctv9yhMxVHjtSNVwEKnn2Q2KdfI724hqI20VTdfg35t94fnoxaIo6QKi4RuQT4ExANPGOMecQr/BbgUWCP4/WEMeaZUMpksViCJ3PoHewdegevz4f+/SErw3e8zonF7D/3XdLPfZc7knbz2OefUrLmfE6uPJ+HHoIHH4SMDLj4YshJeotJCz4g06gpq4ziGioef4UCsMrLEhQhU1wiEg38FbgI2A0sFZFZxpgvvKK+aoy5K1RyWCyW0HN96lX8rXQm5VFCVBQk9VxLeve13PLVE/RuN5jPPoOlS+H99+HYsak8wlT6x3/BuOTPOD95HkMSVxL79GtgFZclCAIqLhH5PvCSMeZYA9MeBWw2xmx10pkBXAl4Ky6LxRLhDO71EN/dAv86oRMy2lXGcHXitYzpfzcAV1+trroadl37VeaXjuOz0vE8d/RW/nH0DpKjShidtIivPwmTJ+t2Qx5WPjmbZUvLKY9pQ0J1MSNGJjD0zinhyKalmSDGmPojiPwamAYsB54D3jeBLtLrvgZcYoy5zTn/BjDa3bpyugp/CxwCNgI/Nsbs8pHW7cDtAJmZmSNmzJgRVOb8UVJSQkpKylml0VREkqzV1VBaWkJsbGB5y8rUMoNni/r6OHFCZ795xkf8YYzGTUsLTtaamhLi4wPLWl6uv762H/GmtFTzFBsbOO7x45CaGng/uNpaqKoKTtbKSs1bUlLg+588CVFREB8fOG5JCSQmBn4GABUV9ctau/MLYmr0FVJ6Mo6VGzuz7IvuFK7rxqGjqQB07lzGqFFHGZy7h149D5OQWFuXgDGkJEFihyAedAAi6f/VGmWdMGHCMmNMnrd/QMUFICICXAx8E8gDXgOeNcZsqeeaYBRXO6DEGFMhIncA1xpjJtYnS15eniksLAwoc30UFBSQn59/Vmk0FZEk66FDsHBhAZ075weMu2KFjnkEs+P03LkwdKhac6iPykooKNBxlEAcOAClpQX07BlY1g0bVCn26xc43SUBdkB2U99+XG5KSnCM7AaWddeu+i1nuAk0Hd7NfGeMK8PPGJebTZsKyM31L+uWDx6kz8uvEF9d51ceHcX6q79JTZ//Zf58WLZM5Ssvh+ioWnpkFtOvy1H6dTlKt/YnSKop4tZnvxpYmABE0v+rNcoqIj4VV1BjXMYYIyL7gf1ANdAWeENEPjTG/MzPZXuArq7zLtRNwvCke8R1+gzw+2DksVgskUuvi+9nI9Bp5utkHK/maHI8ey+/kT5X/C+gCvK227RV/sbvl7NudzvW7c7gnaW9eGdpL5Liq+iTfYzqp3Srlp49Nd2Pp7/Fqg9fo7b6OFExaQy56Bom3jI1bPm0hI5gxrh+CNwEHEaVy0+NMVUiEgVsAvwprqVAroj0QBXWNOB6r7Q7GWP2OadXAOvOKBcWiyWi6HXx/RzNu58562DsWOjtI05SEpzTYTt9Ohdx5egtnDgZy4Y9bVm/O4P1u9py550aLycHBnTfTkbFJnp3hJQEqK0+zop3pwNY5dUCCabFlQF81Rizw+1pjKkVkcv9XWSMqRaRu4D30enwzxlj1orIg0ChMWYW8AMRuQJtxR0FbjnDfFgslhZIt5wEtuyqpDY6jtTEKvJ6HySvxx6ys+MozZ3AwoW6A/THCzpQXvVTALLTt9K742p6d1xF5ey3reJqgQRUXMaYX9YTVm8LyRgzG5jt5Xe/6/ge4J7AYlosltZIu8umwH9ns33bcariUomuLKNT9xS6fGUCoN2K3/oWzHnkOvYc682mA0PYfHAICzZP4dONU4mSGt4cCZMmwYQJ2rpLToZ1n33CZzP+SZsRY3n69Rc4f9pN9D9/QphzawkWaznDYrE0a9pdNoWiLTrmdc45vuPExKbQrd1GurXbyIUDXqeqJpYdR/qx8cBo9pdN5dFH4ZFHdCbroL7FdJCj5GR0YezAGE4cPsQHTz8BYJVXhGAVl8ViiXiyh13DnsLp6KgDxEZX0avjes696DwGTYbiYp3xuWIFLPnkEKsPX0mtuZrn5tfSuc259Oywlk1HV/HQoAm0bRvWrFiCwCoui8US8eROnArA3hWvUVtznFrSyRxyDYMmXwHoUoqLLlJXwI+pqIpn+5H+HEq9mpWLY/ls0xUUbIjlb+20VXf++erGjYPOncOYMYtPrOKyWCwtgtyJU8mdOJXly6FDB+ja1Xe8+LT2cPwQfbNWcOHk7oxLfouq6jgO1YwifcTdfPopTJ8Of/2rxs/J0bGxMWP0d9Cg4BZiW0KHVVwWi6VV0XP8TWx47wlqqyu+9EtMEm67fTT9z9fz6uq6bVs++wzmzIGXX9aw1FQYPRrOO0+V2ejR2O7FJsYqLovF0qrIHKATMLZ++k8AUtt14PzrTp1VGBNTt23LD3+ollO2b1cLIgsWwMKF8PDDao4LoG9fOPdcVWKjRsHgwcGZ/bKcGVZxWSyWVkfmgAlkDpjAnj0FTH3stoD2GkXU8G+PHnDjjep34oRavF+8GBYtgnffhRde0LCEBBg2TPcg87jcXLUNaTl7rOKyWCyWMyA1FSZOVAfaKtuxQ2cvLl6sv888A3/+s4anpaky87TkRoywyuxMsYrLYrFYGgERnciRkwPXXKN+1dWwbh0UFmrrbNkynfRR4QyvpaTAkCGq0IYNU2PSAwcGZ7G/NWMVl8VisYSImBidXn/OOfDNb6pfVRV88YWaqvK455+HJ56ou6ZfP1VoQ4boeNngwdqisyhWcVksFksTEhtbp5Q8yqy2FjZvhpUr1a1apVv0eGYyAqSljWX4cJ2OP3BgnQtmq5mWhlVcFovFEmaioqBPH3WebkaAI0d0X7I1a+D99w9x5Eg2L7ygE0M8ZGXBgAHq+vfX1lr//uofaJPSSMUqLovFYmmmtGsH+fnqzjlnI/n52RijG4Z+/rl2OXqct0JLS9Np+h7Xp49OBsnN1bG1SMYqLovFYokgRHSX7W7dYMqUOn9jYO9eWL9eJ4Rs2KDu00/hpZdOTSMrC3r3ViXWu7duxtmrl/5mZDT/lppVXBaLxdICEFG7ip07w4UXnhpWVqZjaJs2wcaNerx5M7z/vk4McZOWpgrMs27NM1Oye3d1bdo0VY78YxWXxWKxtHCSkupmJ3pTWgrbtsGWLeq2bYOtW7Xl9u67UF5+avw2bepafF27nuq6dGkao8RWcVksFksrJjlZZyoOGnR6mDFw6JCau9q+HXbuVLdjh/4uWqQTSLxJSxvLwoU6YSQUWMVlsVgsFp+IQMeO6kaN8h2ntBR2765ze/bA4sUHycoKXdPLKi6LxWKxnDHJyXUzFz0UFGwiIyN0istaybJYLBZLRBFSxSUil4jIBhHZLCJ3+wiPF5FXnfDFIpITSnksFovFEvmETHGJSDTwV+BSYABwnYh4D9XdChwzxvQGHgN+Fyp5LBaLxdIyCOUY1yhgszFmK4CIzACuBL5wxbkSeMA5fgN4QkTEGGtOMpKprg4cxxi1zxZMXICamsBxa2qCv79nA8DGltWY4GT1UF3d+PkyJvi4jf0MPAQTryH5Mib4fEHD8mWJPCRUOkJEvgZcYoy5zTn/BjDaGHOXK87nTpzdzvkWJ85hr7RuB24HyMzMHDFjxoyzkq2kpISUCLF5Ekmy1tSovFFRgeWtrlb7bMHsRVRVFfxuspWVEBcXOJ4qoxKiowPL6nnBRkcHTre6WuMFY3mgIfmqqQlOVo8yignik7SmRuUM9hnExASXr2Bl9SiiYMrAo5Ab+xnU1pbQpk1k/L8i6V3QWLJOmDBhmTEmz9s/ImYVGmOeBp4GyMvLM/n5+WeVXkFBAWebRlMRSbJCZMlrZQ0NVtbQYGWtI5STM/YAXV3nXRw/n3FEJAZoA/hYzmaxWCwWixLKFtdSIFdEeqAKahpwvVecWcDNwELga8DHgca3li1bdlhEdpylbO2BwwFjNQ8iSVaILHmtrKHByhoaWqOs3X15hkxxGWOqReQu4H0gGnjOGLNWRB4ECo0xs4BngRdFZDNwFFVugdLtcLayiUihr37T5kgkyQqRJa+VNTRYWUODlbWOkI5xGWNmA7O9/O53HZcDXw+lDBaLxWJpWVjLGRaLxWKJKFqr4no63AI0gEiSFSJLXitraLCyhgYrq0PI1nFZLBaLxRIKWmuLy2KxWCwRilVcFovFYokoWqziEpEMEflQRDY5v239xKsRkZWOm+Xy7+FYrN/sWLAPwpBQ6GQVkaEislBE1orIahG51hU2XUS2ufIxNAQynrGlfxG5x/HfICKTG1u2M5D1JyLyhVOOc0SkuyvMZ30Io6y3iMghl0y3ucJudurMJhG5OdSyBinvYy5ZN4pIkSusycpWRJ4TkYOOWTlf4SIif3bysVpEhrvCmrRcg5D1BkfGNSKyQESGuMK2O/4rRaSwGciaLyLFrud8vyus3rrTIIwxLdIBvwfudo7vBn7nJ16JH//XgGnO8ZPAd8IpK9AHyHWOs4F9QLpzPh34Wgjliwa2AD2BOGAVMMArzneBJ53jacCrzvEAJ3480MNJJzrMsk4Akpzj73hkra8+hFHWW4AnfFybAWx1fts6x23DLa9X/O+j6zfDUbbjgeHA537CpwDvAgKcCywOY7kGknWMRwZ0t43FrrDtQPtmVK75wDtnW3cCuRbb4kItz7/gHL8ATA32QhERYCJqsb7B158BAWU1xmw0xmxyjvcCB4GzXowdJF9a+jfGVAIeS/9u3Hl4A7jQKccrgRnGmApjzDZgs5Ne2GQ1xnxijClzTheh5sjCQTDl6o/JwIfGmKPGmGPAh8AlIZLTQ0PlvQ54JcQy+cQY8ylq1MAfVwL/NMoiIF1EOhGGcg0kqzFmgSMLhLe+BlOu/jibun4aLVlxZRpj9jnH+4FMP/ESRKRQRBaJyFTHrx1QZIzxbIywGwjdPtTBywqAiIxCv1q2uLwfdroTHhOR+EaWrzOwy3Xuqzy+jOOUWzFajsFc25g09H63ol/eHnzVh1ARrKxXO8/2DRHx2P9s6nJt0D2d7tcewMcu76Ys20D4y0s4yrUheNdXA3wgIstEd9FoDpwnIqtE5F0RGej4NWq5RoR1eH+IyEdAlo+ge90nxhgjIv7m/Xc3xuwRkZ7AxyKyBn3pNiqNJCvOV+GLwM3GGGdXKe5BFV4cun7i58CDjSF3S0ZEbgTygAtc3qfVB2PMFt8pNAlvA68YYypE5A60VTsxjPIEyzTgDWNMjcuvuZVtRCEiE1DFNc7lPc4p047AhyKy3mkVhYvl6HMuEZEpwFtAbmPfJKJbXMaYScaYQT7cf4ADzkve87I/6CeNPc7vVqAAGIZaqE8XtVgPvi3bN7msIpIG/Be41+ne8KS9z+nyqACep/G74s7G0n8w1zYmQd1PRCahHw1XOOUG+K0PYZPVGHPEJd8zwIhgrw0BDbnnNLy6CZu4bAPhLy/hKNeAiMhg9PlfaYz5cgcNV5keBGYS2m74gBhjjhtjSpzj2UCsiLSnscv1TAfHmrsDHuXUCQ+/9xGnLRDvHLcHNuEMGAKvc+rkjO+GWdY4YA7wIx9hnZxfAR4HHmlk+WLQQeoe1A2sDvSK8z1OnZzxmnM8kFMnZ2wltJMzgpF1GNrNmhtsfQijrJ1cx1cBi5zjDGCbI3Nb5zgjVLIGK68Trx86aUDCVbbOfXLwP4ngMk6dnLEkXOUahKzd0LHhMV7+yUCq63gBujFvOGXN8jx3VInudMo4qLoTtAyhzmS4HDq+Msf5g3zkqXxo19AzzvEYYI1TiGuAW13X9wSWOBXmdc+fLoyy3ghUAStdbqgT9rEj/+fAS0BKCGScAmxEX/j3On4Poi0WgASnnDY75dbTde29znUbgEub4NkHkvUj4ICrHGcFqg9hlPW3wFpHpk+Afq5rv+WU92bgm6GWNRh5nfMH8Pp4auqyRVt7+5z/zG60i+1O4E4nXIC/OvlYA+SFq1yDkPUZ4JirvhY6/j2d8lzl1JF7m4Gsd7nq6yJcytZX3TlTZ00+WSwWiyWiiOgxLovFYrG0PqzislgsFktEYRWXxWKxWCIKq7gsFovFElFYxWWxWCyWiMIqLovFYrFEFFZxWSwWiyWisIrLYmmmiMhIx7hugogki+7FNijcclks4cYuQLZYmjEi8mvUKkkisNsY89swi2SxhB2ruCyWZozozttLgXLUfE5NgEsslhaP7Sq0WJo37YAUIBVteVksrR7b4rJYmjEiMgvdLbYHain+rjCLZLGEnYjeSNJiacmIyE1AlTHmXyISDSwQkYnGmI8DXWuxtGRsi8tisVgsEYUd47JYLBZLRGEVl8VisVgiCqu4LBaLxRJRWMVlsVgslojCKi6LxWKxRBRWcVksFoslorCKy2KxWCwRxf8Hv9jzAW/NKaIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(m.size,1)\n", "fig.tight_layout()\n", "fig.subplots_adjust(top=0.99)\n", "T = np.zeros((m.size,))\n", "absolute_error = np.zeros((m.size,))\n", "relative_error = np.zeros((m.size,))\n", "for j,k in enumerate(m):\n", " nsteps = np.float(k)\n", " dx = 1.0 / nsteps\n", " x_approx = np.linspace(0.0, 1.0, num=k)\n", " T[j] = spi.trapz(f(x_approx), x_approx)\n", " for i in np.arange(0, nsteps):\n", " x_start = i*dx\n", " x_stop = (i+1)*dx \n", " y_start = np.exp(-x_start**2) \n", " y_stop = np.exp(-x_stop**2)\n", " ax[j].fill_between([x_start,x_stop], [y_start,y_stop], facecolor='b', edgecolor='b', alpha=0.2)\n", " ax[j].scatter([x_start,x_stop], [y_start,y_stop])\n", "\n", " ax[j].plot(x, y, 'b-')\n", " ax[j].set_title(r'Trapezoid Rule, $N$ = {}, Approx: {}'.format(k,T[j]))\n", " ax[j].grid(True)\n", " if (j==m.size-1):\n", " ax[j].set_xlabel('x')\n", " ax[j].set_ylabel('y')\n", " absolute_error[j] = np.abs( np.sqrt(np.pi)/2.0 - T[j])\n", " relative_error[j] = np.abs( (np.sqrt(np.pi)/2.0 - T[j])/(np.sqrt(np.pi)/2.0 ) )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tabulate the data and display it" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Results for $e^{(-x^2)}$
OrderRelative ErrorAbsolute Error
40.1650150.146240
80.1587120.140655
160.1576070.139675
\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "results = np.vstack([m.T, relative_error, absolute_error])\n", "import pandas as pd\n", "df = pd.DataFrame(results.T, columns=[\"Order\", \"Relative Error\", \"Absolute Error\"])\n", "df.Order = df.Order.astype(int)\n", "display(df.style.hide_index().set_caption(\"Results for $e^{(-x^2)}$\"))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "jupytext": { "cell_metadata_filter": "-all", "main_language": "python", "notebook_metadata_filter": "-all", "text_representation": { "extension": ".py", "format_name": "light" } }, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }