{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Import libraries\n", "import time\n", "import numpy as np\n", "from numpy import linalg as LA\n", "\n", "import matplotlib.pyplot as plt\n", "plt.rcParams['text.usetex'] = True" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## MA934\n", "\n", "## Finite Difference Method primer\n", "\n", "Here we expand on the content of Lecture 3 and build some of the functionality required to compute numerical derivatives." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "We consider a well-known setup, with\n", "\n", "$$f(x) = \\cos(x),\\ x \\in [0,2\\pi]\\ \\textrm{ (and hence } f'(x) = -\\sin(x) \\textrm{)}.$$\n", "\n", "You are encouraged to consider other tests (e.g. polynomials) in order to explore the implementation." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAv8AAAHzCAYAAAC30XLEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAACMBklEQVR4nOzdeVxVdf7H8dcBF1JUXHChbBKsaZmpBCvbMaF9fjUNLu2rYLbrCLmBoGhg2mJZ4NQ0rSJaVjOVgUZjM1kq1TStBk6WpoRKeptQwe/vD+5lAC/Ies9d3s/H4z6ye8695305LJ977vf7+VrGGERERERExP8F2R1AREREREQ8Q8W/iIiIiEiAUPEvIiIiIhIgVPyLiIiIiAQIFf8iIiIiIgFCxb+IiIiISIDoZHeAQNKvXz9zzDHHePy4P//8M927d/f4ccU9nQ/vovPhXXQ+vIfOhXfR+fAu3n4+Nm7cWG6MCXe3TcW/Bx1zzDFs2LDB48ctKioiNjbW48cV93Q+vIvOh3fR+fAeOhfeRefDu3j7+bAs69vGtmnYj4iIiIhIgFDxLyIiIiISIFT8i4iIiIgECBX/IiIiIiIBQsW/iIiIiEiAUPEvIiIiIhIgVPyLiIiIiAQIFf8iIiIiIgFCxb+IiIiISIDQCr9ebN++fezatYu9e/dSXV3d6ufp1asXX3zxRTsmk7bQ+fAuOh/exVfOR3BwMD169KBPnz507drV7jgiIs2m4t9L7du3jy1bttC7d2+OOeYYOnfujGVZrXquvXv30qNHj3ZOKK2l8+FddD68iy+cD2MMBw4cYM+ePWzZsoWjjz5abwBExGf4VfFvWVYCMNYYM7qFjwsDpgIlzruijDEprd2vPezatYvevXvTr1+/jnh6ERFpJcuy6NKlS+3v5127djFo0CCbU4mINI9fFP+WZeU4/xkJ9GnFU+QDScaYUufzRVqWVWCMiW/lfm22d+9ejjnmmPZ+WhERaUc9e/bkP//5j4p/EfEZfjHh1xiTZIxJoqY4bxHnpwWlroLe+Xyldba1aL/2Ul1dTefOndv7aUVEpB117ty5TXOyREQ8zS+K/zYaC2x0c38BkNSK/dpNa8f4i4iIZ+j3tG9zOBykpWcQHjGYoOBgwiMGk5aegcPhaNV+Ir7AL4b9tFEckOPm/lJgeCv2ExERES/ncDgYcW4sJQd6UTlqJvQbQnn5ZrKXLWXFytdYt7aI0NDQZu8n4isC+sq/cwJvGLDLzeYK57Zm7yciIiL2a86V+vkLFtYU9BdOg/5DISgY+g+l8qLplOzvyfwFC93vB273a+5xRewW6Ff+Dzs52Fn4N2s/Y0xFO2QSERGRVmrulfrFOUtqttcdumUMVO6h8lcjePiRR+nVI5Ts7PlU9omC5xJh7w7Y9zN06QYhPans1IV5895n01df0rNnT15e+Rq7rVCqoi6Cs0+gvHMI2ctW6BMC8SqBXvx3OMuyEoFEgAEDBlBUVNSsx/Xq1Yu9e/e2S4bq6up2ey67VVRUcPTRRze5z8iRI3n11Vc9lKjl2nI+rrjiCjZv3sy//vWvVh/f3XM8/PDDpKam8uqrrzJy5Mh6+69cuZKFCxfy8ccfExYWxr/+9S/CwsJITU1l5cqV/Oc//+HKK6/k2WefbXUmO/nTz0d7a+r7oqM053zYkasplZWVzf7d7kscDodPvq4//+VZNu3rwf6Lp/2vsHdeqd/05mzuuOtubr7xBnZu3wZ9joYtH8E378F3n8CeH+BAJQB7gMmTJ9c8vtse6DkQjjoZjugF+xzwy0/w3584sGsL7777LuXl5ezfv79m/+2bav7bqSuVR53Ml+UVXHPd9dx3z92tnifiq+fDX/ny+VDx38GMMblALsDw4cNNbGxssx73xRdftNtCN76waE5zubpqREZGkpJSs8RCRUVFvX0iIyO9+vW25XwEBwcTFBTUptfn7jkGDBhAZGQk3bp1q3f/8uXLueGGG0hMTGT69OmUlpYyePBgRo8eTWFhIVOnTiUyMhLAq7/mTfGnn4/21tj3RUdqzvmwI1dTQkJCGDZsmN0x2l1RURHN/ZvlTUZfcz37G17RB7As9g+/mjdWzWFMwh/oEhLCvpyEmiv5wV1g8Klw9DDoNRDMQcI+f4XSr77g2BN/y864+/835Keusk2EH9jO1q1bCI8YTPnIadBzAFTugd3fw7fF8O0Gqnd/x+uvlrDhww+Ij4/nvPPO4/Mvv+LZF15i545t9B0QwcSk8UyZPKnRTwd89Xz4K18+H4Fe/Lsbw1+PMaaiOe/SNeTHsyIjI0lMTLQ7ht9ITEx0+/WcN28eCQkJ5OT8b657aWkpy5cvJz8/n4SEdu9yK16kse8LT8nNzSUyMpK4uLh699udS7zbzh3boN+Q+nfu/wX+8yF8/XfKf/ieyy+/nK5duxIUEsbB+D/CkNOg8xE1+xpDyKpM7r5jIr179+aOCYlkL1tK5UXTDxkiFFKcx+2Jt/3vuK65A93Caj5ViDqrZt/d27CeuYHzzjuPv/71r//7pLT3URAzlvLBp5C9bJWGB4lHBPSEX2fBXkHN4mANRTq3NXs/EX9TXFxMnz71p7yUltYsdeG64i/SUVJSUsjPb/HyLeLHmjOhtu+ACCjfXDN+/9uN8OpMeOL38NcM+HYjId268+abb7Jt2zZOGHIUIaV/r7lKX10FZZsIWZVJVJc9TJk8CYApkycR1WUPIasyoWxTo/vVHtedAz/Tb9BRLF26lDvuvIsukcPh7FshtB9syIOVM6i0urJpb3C9CcQiHSGgi3+nDbif0BsFFLZiP5+kDgUiIuLNXBN5s/PfpXzUTMzdb1E+aibZy4oYcW5s7d+riUnj6bL2ScifBCuSYftXcPLvYPRCukbGkJyczMUXX0yfPn1Yt7aI5DGxhK/JJGjRpYSvySR5TGy9q++hoaHN2m9i0nhCipfWvOmoq8EnBE8seYr9Z42HM66B0Qvglmfht5fBl2vYX7KeBx54gHXr1nnmiyoBScV/zarA8W7ujwPyWrGfz2nuL1RfVFFRQVJSElFRUfTu3ZvRo0fXXrl2Wb58Ob1796a0tJTCwkLi4+NJSkoiKSkJy7IOmVNgWRYxMTH17istLcWyLLKzs6moqCAlJYWoqCgsyyIqKqp2fkJzjutSWFhITEwMvXv3JiYmhuXLl7Nr12FHqtXT3OfIzc3FsiyKi4sByM7Orp2U5trWu3fveo+JiYnBsqx6r62iooLRo0fTu3dvt6/7cK+5uY8vLi6u/Rq7zmtjXOe/7rmoe04Pd8ymNPdcu85DcXEx2dnZtV+7+Pj4Q76/mrNvc76OTX3fp6SkYFkWubm59Y5dXFyMZVm1X8+G3xd1j11cXExSUhK9e/emd+/etccvLS0lPj6+9uvR8BjN+ZqNHj269mfPlcH189VYrua87ob5m/s9JN6hOa05P/jgA977+7vs//Yj2P4lDB8LNz0DJ11IyGevM7Tr3tor9VBT2KenpVK2dQvV1VWUbd1CelrqIcNumrNfcz8hOGRYUq9BMOpuGP8inH41+yt/4cwzzyQ2NpY333wT0/DNhEgb+VvxH0YTPfctyyqxLKveKr3OCbmRlmVF1tkvGthljFne0v18UXN7Hfua0tJShgwZQmFhIUlJSWRlZVFaWkpUVBSFhfU/rKmoqCAnJ4f4+Hh27dpFfHx8bTGwbNmy2v2Ki4sJCwujuLi4XtG2fHnNt0BcXBzLli2jsLCQhISE2nHx2dnZ9Yqzpo7rOk58fDzFxcUkJiaSlJTEvHnzDil2mtKW50hMTKSgoKD2NRUUFLB69ep6+2RlZVFQUFCv6BsyZAjFxcVkZWWRlJTk9nU39pqb8/hdu3bVFuuuInL48OEsX778kAKyoqKCqKgoli1bRlxcHDk5OSQkJJCbm8uGDRtqj3nyyScfNnNjmnuuKyoqKC4uJiYmhvXr1zN27Fiio6NrC/3W7Hu4r2NT3/dZWVnExcWRlJRUrzgeNWoUkZGRTQ61cZ2DUaNG1T7X8OHDyc3NZfTo0cTHxxMfH09WVha7du0iKSmp3vdcc75mru8t+N/3X0FBQZPj/Jv7896S7yHxLotzllA5bKzbibyVx5zJA/PmMWLECP71r3+RmZnJ1OQphG/7J0FPXOn2Sn17a+4nBI0OD+rWG447j74Dj2TBggV88803XHrppQwbNowXX3yxtuGFSJsZY3z+BmRRs/rubsBQc5U+B4husN9GoMDN48Ocz+Fqy5nVyHGatV9jt5iYGNNcn3/+ebP3PZw9e/Y0ub3foKMM1+UYJq0+9HbdkyY8YnC7ZWmr3bt3G+c5dnsLCwur3TcuLq7e/7tER0ebyMjI2v/Pz8+vfXx+fn69fQGTkJBQ+//JyckmISHhkH0bO1bd7TU/bv87H00dNzo62gCmpKSk3v2RkZH1sjelJc+Rk5NjALNx48Z69zd8/cYYU1BQYABTUFDg9jXu3r37kOd13dfUa27O413/HxcXd0jOhq8pMTHR7etv6TFbqu65dnG97pycHLf71v1aNmffw30dm/N9v3v3bhMWFlZ7n+v7uu7Xy933heu+ht8XYWFhh+Rxfa9kZWW5/2I1eG3G1P99BZjExMRD9neXq7mvuyXfQ83Vnr+vvck777xjd4R6rKAgw71v1/8bdcNThmPPrf15yMzMPOzfPLulzko3ISeONNxXWP+13FdoQk4caVJnpRtjjNm3b5954oknTL9+/WpfX88+/UzqrHSzd+9em1+FeNvPR0PABtNIPeoXV/6NMSnGmCRjTG9jjGWMGe38/+IG+8UYYw4ZumOMqXA+R67z5vbyT3P38zVuOyO49B1S0wvZy7iuTrq7AbXDIaZOnXrIY11XBF1X610SEhIO6V4TFxdX76phYWEhp512GtHR0bVXJl33N+xIUlfdq7INNTyu68pvYmLiIZNqw8LCGj1GXe3xHC1RUVFRezV3165dlJaWUlpaWnvshp+0uHvNLXm86+vpEh0dfchwptzcXBISEhqdmOw65pVXXtmsYzZXU+e6YRZXFyV3V9qbs2/Dr2NLvu/DwsJYvXo1paWltUPCcnJymj2Re+zYsW7z1s0zfPhwAHbu3NnkczX1NWuO1vy8N+d7SLxLvSvme3+EN+bCs7fVTOr97WX0HXgk06ZN84r2r01p7vCg/fv381jOn9jb/7dw/u3QK4I9u8rJfCSH084616eH5Iq9Ar3Vp1DzC7W8fLP7HsY7N9N3YITnQx1GZGRkk20mXcMMoqOjD9nmKkjWr19f7zkaFjNQUyAUFhZSXFxMdHQ0xcXFLFmyhJKSEpYtW0ZOTk5tkVj38a5io6CgoLaobEzD47r2jYqKavQxh9Mez9ESrmE0y5cvP6TIqpvHpeFrbunj3Z3Xulznv2GB5+6YK1euZOXKlYc9ZmNacq4bchXMzXmMu30bfh1b+n0fHR1NTk4OSUlJJCQktKh9ZsM3kX369Gn2G8u2fM3cac3P++G+h8SzHA4H8xcsZHHOkkb73k9MGl/TcvOYs2H1I1C9H04bCzGjCfn7Y9wxwTfav7qGB81fsJAncjPZuX0bfQdGcHvibfVer2tI7r6LnIuVnXw5/PPPVG9cwZd7d5I04XZeeP45e1+M+CS/uPIvbdPcDgW+qOFkSqDRq3vurni6ioXCwsLaIj86Orp2vHBpaWm9sclQU8BGRUWRl5fH6NGjyc/PJzk5udGMjV1pbY8r9B1xlb8pBQUFbj9ibPj6G3vNzX18w/ajbfHqq68265jutPRct7emPtVoqLHve9f3b0vmk0Drz0FHfs1a8rrb83tI2qa5TScSb7uVrj98Am/MgR794Jon4NexhPz9sXpXzH1BcyYQHzLHoXNIzScA4x6F7r158YXnufnmm9m9e7dNr0J8lYp/afZHkL7EdVVv/fr1h2xzFTmnnXbaYZ8nMjKSsLAw1q9fT0FBQW2B7/qv601BdHR0baE9fvx4EhIS2LhxI4mJiURHR9O3b99mZ3cVdBs3bjzMnh37HC3hurra0gKyvR7fkOv8N/X6Xcf85JNPWn2ctp5r1xXv5lyFbs6+Lf2+z83NrV2wbdeuXR7peNPSr1lzhuK018+72KM5TSfWrl3L2WefjWPPT5wfG0u/zgcIen68Ryby2qXRIbkRJ8J1SwCL5557jhNPPJFXX33V4/nEd6n4l2Z3KPAlkZGRREdHk5ube8jVwJSUFMLCwpq9Ou2YMWNqi/y6w0hc4/6Li4tr3wxUVFRQUVFxyBVZd0VJY8LCwtxmd33S4KnnaImwsDDi4uKYN2/eIV9v19ekIx/vjquzT8PX6/oauI65cOHCVh2zNee6YRZXhxt3Q85asq9LS77vS0tLSUpKIjExsbbzzvLlyw9pzdmeWvo1CwsLa9a5b8+fd/G8Jrv4nPIH5j/4IOeffz6dOnXivffeo+idd/hx23dNtub0B00uGvbTVsIjjmL9+vUMGDCAK6+8kquvvpoff/zRsyHFJ6n4F6D5vY59iWti5JAhQ8jOziY7O5uoqChKS0tbtGqoq796aWlpvQKi7mRgV0EWFhZGZGRkbevC5cuXM3r0aLfj2JuSlZVVmz03N5fc3FxiYmJaVAS3x3O0hGtCquvrnZubW9sHvjmTZ9v6+IaWLFlCZGQkMTExpKSk1Hs+1/loyzFbc65drURd58I1ydnd1fyW7FtXc7/v4+PjiYyMrP0axMXF1baEba9PYBpq6dds+PDhFBYW1u7fVAvW9vp5F89r9Ar3zm+h8GF+cezl1ltv5eOPP2bEiBGeD2iT5gzJHTZsGOvXr2f27NmsWLGCE088Ud/vcniNtQHSrf1v3trq05e4Wn02bNPX1P4JCQkmLCzMhIWFmbi4uENaP7raJjZsc9nwmA3bCG7cuLHR+11tNiMjI01iYqLJysoykZGRpqSk5JBWn40dt6CgwERHR5uwsDCTkJBgCgoKTEJCQrNfe0ueoz1afRpT87VKTEw0kZGRBjDR0dH1WlY252vd1OMby+l6jc15voZtJ7ds2dLkMZtyuHPd8HXn5+ebxMTE2u/H5OTkQ56zOfs25+vY1Pe9q62nu8dHRkaasLAws3v37iZbfTZ8rLtWm66fnbrZm/vzYUzN91rd1+D6nmssQ3N+3lv6PdQcavXZNoe0m76v0HDBXYbgLoauoaZn734eyeFt9u7da046NaamLeh1TxruWWW47kkTcuJIc9KpMfXafe7du9dMuH2i6dS5swHMEd1DzYzUNLUE7UC+3OrT9oI4kG4q/sUYnQ9v44nz4SrW3b1hasu+/sgXfz5U/LdNvb73E1YYjjmtpq/9MaebrsedVdv3PhDt3bvXpM5KN+ERg40VFGTCIwYf0ue/9k3CSRcYrllsOPlyA5ig0H7m+N+eqjcAHcSXi38N+xERERHbuJpOdHl9Ojw/Ab77BE4bR9cjujG02z6fbDrRXuoOyV2zerXbIbn1JkwP/DWMuhdG3snBn3fx1dffMGNmqn0vQLySin8RERHpEA6Hg7T0DMIjBhMUHEx4xGDS0jPqLVAVGhrKooXz6bTtU6z/7sQ6WEX41n+QMnakzzad8KRDJkxbFgz7PVyZiTlYxaOPPsIHH3xgb0jxKlrkS0RERNqdq39/yYFeVI6aCf2GUF6+mexlS1mx8rXawv61115j3LhxDBo0iLfeeotjjz3W7ug+pdEJ00NOh7GPYp5P5Pzzz+eZZ55h3Lhxng8oXkdX/kVERKTdNad//5NPPsnvf/97fvOb3/D++++r8G+FJluCcpC+AyI47bTTuPrqq0lLS6uZ8CkBTcW/iEgHS0hIwBhTux5Ee+0r4s2a7N8/bAzzH3yQ22+/nUsvvZR33nmH/v372xPUxx2uJegdtydRWFjIjTfeSEZGBuPGjeOXX36xJ6x4BRX/IiIi0u4aHY5SfQA2LucXx14SExN55ZVX6N69u+cD+gnXhOmQVZlQtgmqq6BsEyGrMonqsocpkyfRtWtX/vznP5OVlUV+fj7nn38+P/zwg93RxSYq/kVERKTduR2Oss8Br0yDL1fTrUcvnnzySTp10vTDtggNDWXd2iKSx8QSviaToEWXEr4mk+QxsfUmTFuWRXJyMi+//DKfffYZI0aMYPPmxoYLiT9T8S8iIiLt7pDhKHt/hLz74PtP6BxxPH+cPAmr4ZAgaZW6LUGrq6vctgR1ufLKK1m7di179+5l5MiRfPvttzYkFjup+BcREZF2V284yqa18NKd8NM2Oh/1G47r3z2g+/fbxdV69aLLr2B3RQVbtmzh1FOH8cUXX9gdTTxIxb+IiIi0O9dwlMT4Uwj622xwlNM7rBdTb7lK/ftt4Gq9mp3/LuWjZsI9qzAX3U/Fnr2ceuowvvrqK7sjioeo+BcREZEWac7iXQB79+7lb399nZ49e/DJJ5+wq2x7o8NRpGO5bb16YhyMeYj9VVWcddbZbNu2ze6Y4gEq/kVERKTZGl5BNne/RfmomWQvK2LEubG1bwDKy8uJi4tj+/btvPXWW5x88sk2Jw9sjbZejTgR4iexe/cuRo4cqS5AAUDFv4iIiDRbcxbv+umnn7j44ospLS3l9ddf54wzzrA7dsBrtPUqwPFxgMXWrVsZNWoUO3bs8Gg28SwV/yIiItJsTS7eFT2WxU/mcvnll/PJJ5+wfPlyRo4caU9QqafJlYB3bqbfoCN54403+Pbbbxk1ahRlZWWeDSgeo+JfREREmq3JK8i9jqR8+1b++c9/8sILL3DZZZd5Npw06nArAd+eeBvnnXcef/vb3ygtLSUuLo7y8nJ7wkqHUvEvIiIizdboFeSD1fDaTACWLFnCmDFjPJxMmtKclYABYmNjef3119m0aRNxcXHs3LnT5uTS3lT8i0+pqKjAsqzD3gJZbm4ulmVRXFzc7s+dnZ2NZVkUFha2+3OLiG9wewXZHIRV2fDdx1x08cXccsst9gUUt5q7EjDAqFGjePXVV/nyyy+56KKLDuniJL5Na2qLT4qMjCQlJQWoeUMgnhEWFkZkZKTdMUTERlMmT2LFytcoWZVJZfRY6HMMvDUXvv47/QdGsDw/3+6I0gjXSsDpaamH3ffCCy9kxYoV/N///R/XXnstL7/8MsHBwR5IKR1Nxb/4pMjISBITE+2O4VZubi6RkZHExcXZHaXVGnsNiYmJXvt1FxHPcF1Bnr9gIYtz5lD+w1bAcNbZZ/PWm2+qh78fueyyy3j44Ye5++67SUlJ4cEHH7Q7krQDDfsRaWcpKSnk+/iVL394DSLScVxXkGek/BEwTJgwgffWrqVHjx52R5N2UHcRt3vuvZeQbt1ZsGABixYtsjuatANd+RcREZEWe/vtt5k0aRJXXnkljz/+eMDPt/IXrkXcSg70onLUTOg3hMqybwh6dQZ33303xxxzDL/73e/sjiltoCv/4rdSUlKwLIvc3Nx69xcXF2NZFqNHjwZq5gykpKQQFRWFZVlERUXVzidwJykp6ZB9KyoqGD16NJZlUVFRUTvp1rIssrOzm5W3uLiY+Pj42sdFRUUd8tiKiora4/fu3ZvRo0dTWlp62OdubBJwTEwMvXv3rv3/w72Gxp6nObmWL19O7969KS4urv16u/YVEd+yadMmxo4dy0knncRzzz1HUJDKCX/hdhG3gb/m4I1/xurandGjx/DFF1/YHVPaQFf+fdS9997Lxx9/3Kx9q6urvWKSzqmnnsrDDz/cLs9VWlraaMeZyMhIIiMjycrKori4mKSkJOLi4monqo4aNYrIyMjaYS3Lli2jsLCQhIQETjvtNNavX092djYVFRXk5OTUPm9FRQUxMTHs2rWLMWPGEBMTQ0lJCbm5ucTHx5OVlUVSUhLx8fHExcXVvoEYPnx4s15T3Vy7du1i48aN5OXlkZycXPuaY2Ji6NOnD0lJSYSFhZGTk0NUVBQFBQXtMsegNa+hubl27dpV+ybJ9dz5+fksX76clJQUsrKy2pxfRNrG4XA4x/IvYef2bfQdGMHEpPFMmTypdiz/Tz/9xP/93/8RHBzMq6++qjH+fmZxzpKaK/4NP8kJCcVcMo39r83k8ssv54MPPqBfv372hJQ2UfEvPqm0tJT4+Hi327KysmoL5vz8fIYMGUJ8fDwlJSWMHj2aiooKNm7cWLt/w0msCQkJFBcXk5ubW6/4T0lJobS0lJKSknodb+oWra77Wzrht7CwkIqKCrKyshp9XFJSEgAlJSX1ssfExJCUlFTv/tZyvXFy/bs5r6GluSIjI2u/romJiViWxfLly1X8i9jM3XCP8vLNZC9byoqVr7FubRFHHHEEV199Nd988w0FBQUMGdLIYl/is5pcxO1Xw8HA1q1b+f3vf09hYSFdu3b1bEBpMxX/PqolV9D37t3rd5Ow4uLiKCgoOOx+YWFhrF69mpiYGGJiYiguLiYnJ+ew7Srj4+MpLCyktLS0dt/c3FwSEhI6pNWl6zld2Roew/VJh7sCOSsri/j4eJYvX05CQkK7Z2tKa3I1fNMWHR3drKFLItKx6g33cF317T+UyoumU7Iqk/kLFvLLzw7efPNNnnjiCWJjY23NKx2j74AIyss31wz5aWjnZvoNOpJFC+czbtw4brvtNp599lnN9/AxGqQHWJalxuV+LDo6mpycHIqLi0lISHDbqrK0tJTs7Gzi4+Pdjvl3jXFv7NOGtnK1Ll2+fHntWPi6GVzHj46OPuSxriE569ev75BsTWlNLnf7ioj9FucsoXLY2EOHe1gWldFjeejhR5g/fz4TJ05kwoQJ9oSUDud2ETcAYwgpzuP2xNsYO3YsGRkZPP/882RmZtoTVFrNb4p/y7LCLMvKsiwr0XlryRiCJMuyCizLSrAsK86yrOiGtzrHKXE+f5jzlmBZVo5lWWHt/6qkvbg+JXC36q2r4M7Ly2P06NHk5+fXDhvypJycHAoKCmqvkmdnZxMTE1NvH3cLmu3atcsT8ZrUklx9+vTp4DQi0hpNDvc4sI+9FbuIjY1tt7lb4p2mTJ5EVJc9hKzKhLJNUF0FZZsIWZVJVJc9TJk8CYAZM2Zw3XXXMXPmTPLy8mxOLS3hN8U/kA/kGGNyjTG5QI5lWYcfF1IjEohzPkcBsLHBLanBvjnAbuctC8gyxlS0x4uQ9pebm8vy5ctrJ9I27C4zfvx4EhIS2LhxI4mJiURHR9O3b996+7iuVtedK9CU1hbkcXFx5Ofns3v3bpKTkykuLqa0tLT2+O6u7rve0Jx22mmtOmZjmvMa7MglIh2j74AIKN986Ia9P8JrqQQFB5Ofn0/nzp09H048xrWIW/KYWMLXZBK06FLC12SSPCaWdWuLaid4W5bFn/70J84++2xuvPFG1q1bZ2tuaT6/KP4ty0oASo0xtQOHXf92bjucUmOM5e4GpBhj6hb/KUAMEA/EGGOi6h5XvEtpaSlJSUkkJiaSkJBQ213G1f6zoqKCioqKQ8bYuytmExISyM3NPWR8ekVFRb37wsLC3F4Jb4orR12uIUaueQfR0dHk5uYesl9KSgphYWHNGu9ft6CvqKhw+0lIS15De+USEfu5He5xYB+8lgb7HCQmJqm7S4BwLeJWtnUL1dVVlG3dQnpa6iGdnbp27corr7zCkUceyZVXXsn27dttSiwt4S8TfsdSc8W+oQJqrtovP8zj3bZJsSwrDjikOjLGuK+YxGNcY/Qb45qYGx8fX6+7TFxcHImJiSQlJTF8+HCio6OJjIysbe0ZHx9PXl4ey5cf+i2zZMkSiouLiYmJITExkaioKDZu3Ehubm69DkPDhw+nsLCQ7Ozs2k43dbsGuVNYWMjo0aNru+RAzYTZsLCw2o47+fn5xMTEMGTIEKZOnVr7vKWlpYed/Ox6c+OamOta26AxLXkNbcklIt5jyuRJrFj5GiWrMqmMHgt9joHX02DHVwweEsX8bHXkkkOFh4fz6quvcvrpp3Pttdfy9ttve0V7cWmCMcbnb9QMv4lzc38CsLsNz5vcnPuae4uJiTHN9fnnnzd738PZs2dPuz2X3Xbv3m2Aw97y8/NNQkKCAczGjRsPeZ7IyEgTFhZmdu/ebTZu3Giio6MNYCIjI01iYqLJysoykZGRpqSk5JDjJyYmmsjISAOY6Ohok5WVVW+fgoICExYWZsLCwkxcXJwpKCiot93d+di9e7fJysqqzeF6rLvjJyQk1Hv+hvvk5OS4fd2ur0fdXImJiSYuLu6QPI29hsaeuy25oqOjTVhY2CEZPMWffj78gS+ej/b8fW23vXv3mtRZ6SY8YrDBsgxgRl4wyuzdu9fuaAHvnXfesTtCk5566ikDmFmzZtkdxSO8/XwAG0wj9ahlGs7m9jHOiba7qRmCU9xgWxxQYGqG77T0eZONMYdcWnbd75wE3McY436lKTeGDx9uNmzY0Kx9v/jiC0444YTmB26CP7b69GU6H95F58O7+OL5aM/f195i1apVXHLJJYwZM4aXXnpJrRy9QFFRkde1V627KFz59q10DTmCfb/8wuuvv8bll19ud7wO5Y3noy7LsjYaY9yu0OkPw34O2zrEsqww04IJuc7CvrFx/Kc55xEUG2OKnV2FSkzNJGMRERGftn37dm644QaOOeYYnn76aRX+4pa7ReH2bf8S6+X7+f3vf8+XX35JVFSU3THFDX8o/jvCVGPM6Ea25dS92m+MSXG2/yx19ymAZVmJQCLAgAEDKCoqalaAXr16sXfv3pYnd6O6urrdnkvaTufDu+h8eBdfPB+VlZXN/t3u7Q4ePEhycjI//fQTqampfPjhh3ZHEieHw+FV32d//suzbNrXg/0X11kULuIkzNhHqHo+iXPOOYfnn3/eb8f/e9v5aAkV/w0cbsGvRob5LKem5WdMww3OTwRyoWbYT3M/Ivriiy/a7aNvX/wY3Z/pfHgXnQ/v4ovnIyQkhGHDhtkdo0l1h2fs3LGNvgMimJg0nimTJ9Xr4JKVlcXGjRvJycnhuOOO8+phDYHG24aZjL7mevaPmnnoonDhkTDiOra//yxr165l1qxZtuTraN52PlrCH1p9HrYZeUuG/FDTHailS6WWAFq2VEREvI5reEZ2/ruUj5qJufstykfNJHtZESPOjcXhcADwwQcfMGPGDEaPHs348eNtTi3erslF4U67BrDIyMhg9erVHs0lh+fzxb+zsK+gZvGthiKd21oigUbG+zuH96hpuYiI+Iz5CxbWjMu+cBr0HwpBwdB/KJUXTadkf0/mL1jITz/9xLhx4zjyyCPJzc3VOH85rEYXhQPY9R/6DYzg+OOP59prr1X/fy/j88W/0wbcT/yNAprdjcfZOaipNwwVuOn77zyOev+LiIjXWZyzhMphYw8dnmFZVEaPZXHOEpKSkvjuu+946aWXCAsLsyWn+Ba3i8IBGENIcR4TJySSn5/Pnj17uOaaa6iurrYnqBzCX4r/fGpW3G0oDshrwfO4Pj1obChRnnG/mm8C0PQqTiIiIjZocnhG3yHs3L6VvLw8Zs+ezZlnnunZcOKzpkyeRFSXPYSsyoSyTVBdBWWbCFmVSVSXPUyZPImTTjqJxYsX88477zB79my7I4uTXxT/zkm1kXUn6zrbde4yxtRbqtU5dGdjI0/lenxFI9sLnd176j5fMlCqVp8iIuKNmhyeUfpPDDBq1KgmV/0WaSg0NJR1a4tIHhNL+JpMghZdSviaTJLHxLJubVHtRPKbbrqJG2+8kYyMDAoLmz0YQzqQP3X7GQVMtSyrxPn/UcYYd58GVND4lf1S583tdmdff5y9/QHCqOnx7+44bWaM0bhLEREv5gsLZU5MGk/2sqVUXjS9/tCfA5VYbz/IEUd047nnniMoyC+uB4oHhYaGkp6WSnpaapP7Pf7443z44Ydce+21fPzxxwwaNMhDCcUdvyn+nRN/D3vZwhhzSDvOOtuKqRm/39Tji/HA+P7g4GAOHDhAly5dOvpQIiLSSgcOHPD6PuZTJk9ixcrXKFmVSWX0WOg7BHZuJvhvs6ne9zPPrVihYkw6VPfu3cnPz+f000/n2muvpbCwUG82baSvvJfq0aMHe/bssTuGiIg0Yc+ePV6/LoG74Rk9/zad6t1bueuuu7jqqqvsjigB4KSTTuLRRx/lnXfeYdGiRXbHCWgq/r1Unz592L17N+Xl5ezfv98nPloWEQkExhj2799PeXk5u3fvpk8fd83mvItreEbZ1i1s3lxKMAcZPnw4Dz74oN3RJIDccsstXHbZZdx///189dVXdscJWH4z7MffdO3alaOPPppdu3bxn//8p00tsiorKwkJCWnHdNIWOh/eRefDu/jK+QgODqZHjx4cffTRdO3a1e44zVZVVcW1117LgQMHeOmllzS0VDzKsiyWLFnCSSedxE033cTatWvp1EmlqKfpK+7FunbtyqBBg9o8FrOoqMjrl54PJDof3kXnw7vofHSs+fPn89577/H8888zdOhQu+NIgHA4HMxfsJDFOUvYuWMboT3DWLduHZmZmaSlpdkdL+Co+BcREQkAn332GbNmzWL06NFce+21dseRAOFwOBhxbmzNKtOjZkK/Iez9sZSg12eSnp7OxRdfzBlnnGF3zICiMf8iIiJ+rqqqiptvvpmePXvy2GOP2R1HAsj8BQtrCv8Lp0H/oRAUDAOO5eA1T2KCOnHllVeyf/9+u2MGFBX/IiIifm7BggWsX7+exx9/nP79+9sdRwLI4pwlVA4bW3+NCYBuYXDObWzfvp05c+bYki1QqfgXERHxUQ6Hg7T0DMIjBhMUHEx4xGDS0jNwOBy1+3z++eekpqaSkJDAmDFjbEwrgWjnjm3Qb4j7jadeCVjMnTuX9evXezJWQNOYfxERER/kbix1eflmspctZcXK11i3toiQkBBuvvlmevToweOPP253ZAlAfQdEUF6+uWbIT0M7N9N3YARHdLK48cYb2bhxI0cccYTnQwYYXfkXERHxQW7HUvcfSuVF0ynZ35P5CxaycOFCPvzwQw33EdtMTBpPSPFSaLhekTGEFOdxx4REnnrqKb744gtmzpxpT8gAo+JfRETEBzU6ltqyqIwey6LHnyA1NZWrrrpKw33ENlMmTyKqyx5CVmVC2SaoroKyTYSsyiSqyx6mTJ7EhRdeyO23387ChQtZu3at3ZH9nop/ERERH9TkWOreR7P7x+2EhoayePFirIZvEEQ8JDQ0lHVri0geE0v4mkyCFl1K+JpMksfEsm5tEaGhoQBkZ2czZMgQbrrppnpzVqT9qfgXERHxQX0HRED5Zvcb3/sTAI899hgDBgzwYCqRQ4WGhpKelkrZ1i1UV1dRtnUL6WmptYW/a59nnnmGzZs3k5ycbGNa/6fiX0RExAc1OpZ657fw8UqOP/4Exo4da084kVY499xzmTRpEk888QRvv/223XH8lop/ERERH+R2LPX2r7CW3UtwkMVf//q6hvuIz5kzZw4nnHACt912G3v37rU7jl9S8S8iIuKD3I2l7v7XqZhf9pCbm0tUVJTdEUVaLCQkhKeeeorvv/+e1NRUu+P4JRX/IiIiPqruWOrPPvs3VfsqueKKK7j55pvtjibSameeeSYTJkzg0UcfZePGjXbH8Tsq/kVERHxcdXU1t9xyC926deOJJ57QcB/xeXPnzqV///4kJiZSVVVldxy/ouJfRETExy1evJj333+fRx99lEGDBtkdR6RVHA4HaekZhEcMpk/fvvxcuZ/i4mIWLFhgdzS/0snuACIiItJ627ZtY/r06cTHx3PttdfaHUekVRwOByPOja1ZtXrUTOg3hL0/lhL06jSmTZvGFVdcwfHHH293TL+gK/8iIiI+bNKkSezfv5/HH39cw33EZ81fsLCm8L9wGvQfCkHBMOBYDo55lIMGrrzy93ZH9Bsq/kVERHzU22+/TV5eHlOnTuXYY4+1O45Iqy3OWULlsLHQ8A1s2CA49Uq++upLXnnlFXvC+RkV/yIiIl6o7vjnoOBgwiMGk5aegcPhAKCyspI77riDY489lpSUFJvTirTNzh3boN8Q9xvPvhWAu+66iz179ngwlX9S8S8iIuJlXOOfs/PfpXzUTMzdb1E+aibZy4oYcW4sDoeDBx54gG+++YbFixcTEhJid2SRNuk7IALKN7vfWPEdYf36s23bNmbOnOnZYH5Ixb+IiIiXcTv+uf9QKi+aTsn+ntw/dRrz5s3j6quvJi4uzu64Im02MWk8IcVLwZj6G4whpDiPu++8g4kTJ7Jo0SLWr19vT0g/oeJfRETEyzQ6/tmyqBw2htzcXI444ggWLlxoT0CRdjZl8iSiuuwhZFUmlG2C6ioo20TIqkyiuuxhyuRJZGZmMnDgQPX+byMV/yIiIl6myfHPO7/lwP59tYWQiD8IDQ1l3doiksfEEr4mk6BFlxK+JpPkMbGsW1tEaGgovXr1YtGiRXz88cc8+uijdkf2WSr+RUREvEyj458rHVC0mE6dOzNhwgTPBxPpQKGhoaSnpVK2dQvV1VWUbd1CeloqoaGhtftcddVVXH755cycOZNvv/3WxrS+S8W/iIiIl2l0/PM/nobKPdx08y0EBwfbE07ERpZl8dhjjwFwxx13YBr+jMhhqfgXERHxMm7HP39eAJ+8Sp9+4Ty04EG7I4rY5le/+hWzZ8/mb3/7Gy+//LLdcXyOin8REREv03D8s/XoJXRavYDQ0FD+9fFH9YZBiASiu+++m1NPPZW77767du0LaR4V/yIiIl6o7vjnRx55mKoDB/jTn/7EkUceaXc0Edt16tSJxYsXs23bNubOnWt3HJ+i4l9ERMSLbdu2jenTpxMfH8+YMWPsjiPiNc4880xuuOEGFixYwDfffGN3HJ/hN8W/ZVlhlmVlWZaV6LxltfDxJc7HhTlvCZZl5ViWFdaexxEREWmJSZMmsX//fh5//HGshn3/RQLcAw88QJcuXZg0aZLdUXxGJ7sDtKN8IMkYUwpgWVakZVkFxpj4Zj4+Eshx3gBKgXhjTEU7H0dERKRZVq9eTV5eHrNmzeLYY4+1O46I1xk0aBCpqakkJyfz5ptvcskll9gdyev5xZV/y7ISgFJXQQ5QpzhPaObTpAAxQDwQY4yJqvt87XgcERGRw6qqquKee+5hyJAhpKSk2B1HxGs4HA7S0jMIjxhMUHAwWQsfoW/fvtx1113s37/f7nhez1+u/I8FCtzcXwAkAcub8yTGmGJPHEdERORwcnJy+Oyzz3j55ZcJCQmxO46IV3A4HIw4N5aSA72oHDUT+g1hZ/lmOr/7ODtL/kV2djYzZsywO6ZX84sr/0AcNcN0GioFhvvgcUREJIDt3LmTmTNncsEFF3DllVfaHUfEa8xfsLCm8L9wGvQfCkHB0H8oBxIWEhTal9mzZ/PDDz/YHdOr+Xzx75yQGwbscrO5wrmtJc8XbVlWXEcfR0REpDGzZs3ip59+4uGHH9YkX5E6FucsoXLYWGj4c2FZHBx1D/v37+f++++3J5yP8PniH+hzuB0aduxpxGnOcfsVxphCV0efDjiOiIgEsIbjlcMjBpOWnlG7UNG///1vnnjiCSZMmMBvf/tbm9OKeJedO7ZBvyHuNx5zBmDx7LPP8v7773s0ly/xlzH/7SHHGFPo+h9jTIqz/Wdp3ftFRERay9145fLyzWQvW8qKla/x/t/f4d5776Vnz55kZGTYHVfE6/QdEEF5+eaaIT8N7dxMv4ERdAky3HXXXXz44YcEBfnDde72ZRlj7M7QJpZlRQIl1HToKW6wLY6aybi93bTsbM5zZwFxxpiY1h7H+elBIsCAAQNili5d2tIYbeZwOLQUvBfR+fAuOh/exd/Px5//8ixLP9jM/otn1B+2YAxd3pzNWRHBFK1Zw1133cVVV11lX1D8/1z4Gp2PGn/+y7MsXVfK/ktmuv0ZGjcikqMiBjF37lz++Mc/ctlll3VIDm8/HyNHjtxojHE7H9Ufiv8wYDdNFOXGmFYNmHQW7jnGGKs9jjN8+HCzYcOG1kRpk6KiImJjYz1+XHFP58O76Hx4F38/H+ERgykfNdP9VcsfPido2b0c/+tf8/HHH9O5c2fPB6zD38+Fr9H5qFH76dn+nlRGj4W+Q2DnZkKK84jqsod1a4vo3r0755xzDps2beLrr78mLCys3XN4+/mwLKvR4t/nPwtxXmmvoGaRroYindua5Bze02Sf/vY4joiIBLYmxytv+YiD1dU89NBDthf+It4qNDSUdWuLSB4TS/iaTIIWXUr4mkySx8Sybm0RoaGhWJbFokWLKC8vJz093e7IXsfni3+nDbifkBsFNGe8fgXgrsd/VIP723ocEREJYH0HRED55kM3OHbCBy/QpWsIF154oeeDifiQ0NBQ0tNSKdu6herqKsq2biE9LbXeMJzo6GgSExNZtGgRn332mY1pvY+/FP/51KzM21AckNeMx+c1XM3XKQHIacfjiIhIAJuYNJ6Q4qXQcMjtP56Cqv2MT0x0/0ARabE5c+bQo0cP7rnnHnx9mHt78ovi3xiTC0Q6J+UCNf36gV3GmHqr7jqH+Gxs8BSFDdp6YllWMlDqfO4WH0dERKShKZMnEdVlDyGrMqFsE1RXwecF8Nkq+vbvzwNzM+2OKOI3+vXrx+zZs1m9ejWvvPKK3XG8hj+1+hwFTLUsq8T5/1HGGHdX6StosFCXMabYsixXdx+oWbCrpJHHN/c4IiIi9bjGK89fsJAncjMp/2ErwZ060bV7d/71UbFXdw8R8UUTJkwgJyeHyZMnc+mllxISEmJ3JNv5TfHvnJCb0oz9Yhq5vxj34/5bdRwRERF3XOOV09NSefHFF7n22mvJefIJIiIi7I4m4nc6derEwoULufDCC3nsscf44x//aHck2/nFsB8RERFf8/PPP5OcnExMTAw33XST3XFE/FZ8fDyXXHIJc+bMYefOnXbHsZ2KfxERERtkZWWxdetWHnnkEa1CKtLB5s+fz969e7VyNir+RUREPO7bb79l/vz5XH311Zx99tl2xxHxeyeddBK33XYbixcvZtOmTXbHsZWKfxEREQ+bPn06AA888IDNSUT8l8PhIC09g/CIwQQFB7N85esEBQUF/Lh/v5nwKyIi4gs2btzICy+8wP3338/RRx9tdxwRv+RwOBhxbiwlB3pROWom9BvCrvLNdFr1AK+99hqrVq3ioosusjumLXTlX0RExEOMMUyZMoV+/fpx//332x1HxG/NX7CwpvC/cBr0HwpBwdB/KFVjF0GnLtxyy60cPHjQ7pi2UPEvIiLiIW+++SbvvPMOqamp9OrVy+44In5rcc4SKoeNBcuqv6HLEXDaOLZt20peXp494Wym4l9ERMQDqqqqSE5OZujQoSQlJdkdR8Sv7dyxDfoNcb9x+DgApk6dSmVlpQdTeQcV/yIiIu2g4eTC8IjBpKVn4HA4AHjmmWf47LPPmDdvHl26dLE5rYh/6zsgAso3u9+4ewu9+obz7bffsmjRIs8G8wIq/kVERNrINbkwO/9dykfNxNz9FuWjZpK9rIgR58ayY8cOUlNTOfPMM/nDH/5gd1wRvzcxaTwhxUvBmPobjCGkOI977rqTyy67jMzMTMrLy+0JaRMV/yIiIm3U2OTCyoumU7K/J2PHXc0PP/zAgw8+iNVwDLKItLspkycR1WUPIasyoWwTVFdB2SZCVmUS1WUPUyZPIjs7OyAX/lLxLyIi0kaNTi60LCpPuJh33y3iqquu4qyzzrInoEiACQ0NZd3aIpLHxBK+JpOgRZcSviaT5DGxrFtbRGhoKCeeeCLjx4/niSee4Ouvv7Y7sseo+BcREWmjJicXfvUuGKMFvUQ8LDQ0lPS0VMq2bqG6uoqyrVtIT0slNDS0dp/09HRCQkJISUmxMalnqfgXERFpo0YnF+78Fj57i5DuoRx77LGeDyYiTRowYAD3338/K1eu5O9//7vdcTxCxb+IiEgbNTq5cO0SsCzumDjRnmAiclj33XcfRx11FJMnTw6Ihb9U/IuIiLSR28mFn74Bpe/Tf8AAZqXOtDuiiDSiW7duZGZmsmHDBpYuXWp3nA6n4l9ERKSNGk4utB69hE5Fi+jRoyeffvxRvTHGIuJ9rrvuOoYNG8a0adPYt2+f3XE6lIp/ERGRdlB3cuGLL75A1YH9PProI/Tv39/uaCJyGEFBQWRlZfHtt9/y5JNP2h2nQ6n4FxERaUf79u1j6tSpnHzyyVx//fV2xxGRZoqPj2fUqFHMmTOHPXv22B2nw6j4FxERaUeLFy/mP//5D/Pnzyc4ONjuOCLSAvPmzaO8vJwFCxbYHaXDqPgXERFpJ7t372b27NlceOGFXHjhhXbHEZEWOu2000hISGDBggXs2LHD7jgdQsW/iIhIO8nOzmb37t1kZWXZHUVEWikzM5PKykoyMzPtjtIhVPyLiIi0gx9++IFHHnmEa665hlNPPdXuOCLSTA6Hg7T0DMIjBhMUHMzZsaM45dRTefLJJyktLbU7XrvrZHcAERERfzB79mwOHDhARkaG3VFEpJkcDgcjzo2l5EAvKkfNhH5DKC/fzN4P/kJVVRXTpk3zu97/uvIvIiLSRiUlJSxZsoTx48cTFRVldxwRaab5CxbWFP4XToP+QyEoGPoPZd/lGQT1GUxe3jI++eQTu2O2KxX/IiIibZSamkrnzp2ZOVMr+Yr4ksU5S6gcNhYsq/4Gy6L6gnuwLJg2bZo94TqIin8REZEmNBwPHB4xmLT0DBwOBwAff/wxL774Ivfeey+DBg2yOa2ItMTOHdug3xD3GyN+A8Abb7zB3//+dw+m6lga8y8iItKIxsYDZy9byoqVr7FubRHTp08nLCyMKVOm2B1XRFqo74AIyss31wz5aWjnZvoOjKBrEKSkpPDPf/4Tq+EnBD5IV/5FREQa0dh44MqLplOyvyd33n0Pb7zxBvfffz+9e/e2O66ItNDEpPGEFC8FY+pvMIaQ4jwmJiUya9Ys1q1bx6uvvmpPyHam4l9ERKQRTY0Hrhw2hhdeeJ5BgwZx11132RNQRNpkyuRJRHXZQ8iqTCjbBNVVULaJkFWZRHXZw5TJk7jpppv49a9/zbRp06iqqrI7cpup+BcREWlEk+OB9/5I1f79pKam0q1bN88GE5F2ERoayrq1RSSPiSV8TSZBiy4lfE0myWNiWbe2iNDQUDp16kRmZiZffPEFzz33nN2R20zFv4iISCP6DoiA8s2HbjAH4e85BAV34tZbb/V8MBFpN6GhoaSnpVK2dQvV1VWUbd1CeloqoaGhtftcddVVnH766aSlpVFZWWlj2rZT8S8iItKIRscDf7EGKrZy5ZVX0rlzZ3vCiYjHWJbFAw88wHfffcfjjz9ud5w2UfEvIiLSCLfjgX/4Amv1Q4QccQR/fvopuyOKiIeMHDmSCy+8kLlz59a2+vVFflP8W5YVZllWlmVZic5bVgsfH21ZVrLzOQosy8qxLCvMzX4lzucPc94SGttXRER8m7vxwKFvzMAcqOSF55+nZ8+edkcUEQ964IEH2LVrF0uXLrU7Sqv5TfEP5AM5xphcY0wukGNZVkFzHmhZVjQw3BiTbYxJMcbEAyXAbsuyIhvsHgnkALudtywgyxhT0V4vREREvEfd8cB79vxE95CunHfeefz+97+3O5qIeNiwYcMYN24cK1asYMeOHXbHaRW/KP4ty0oASo0xpa77XP92bjucOOcbhlrGmGygmJpCv64UIAaIB2KMMVF1jysiIv7r0UcfZceOHcybN88vFvsRkZZLT09n//79zJs3z+4oreIXxT8wFtjo5v4CIKkZj0+yLCvOzf2FwCH3G2OKjTGFxpjilsUUERFftWvXLrKysvjd737HWWedZXccEbHJcccdx0UXXcQTTzzBd999Z3ecFvOX4j8OcHf1vRQY3sznaDi8R0REpFZWVhZ79uwhMzPT7igiYrMbbrgBYwxz5syxO0qL+Xzx75xoGwbscrO5wrmtSc6hO7luNkVSM/TH3XGjG/m0QERE/My2bdt49NFHue666/jtb39rdxwRsdnAgQNJTEzk6aefpqSkxO44LeLzxT/Q53A7tKYTj/MxCUDDAV2nOecRVBhjCl0dhlr6/CIi4jvmzp1LVVUVs2bNsjuKiHiJ6dOn07lzZ5/7vWCZhguX+BhnN54SaibfFjfYFkfNuP/eLe3GY1lWDtDHGDO64XMaYwob3FcCJDW837ktEUgEGDBgQIwdraEcDke9VerEXjof3kXnw7t44/nYvn07119/PZdeein33Xef3XE8xhvPRSDT+bDfL7/8wtJl+bz61zfYs/NHevYN56iB4Xz+2Wc89dRTDBkyxO6ItUaOHLnRGON26HsnT4fxBc43DcONMTENt7kr8IHl1LT8dLd/LpALMHz4cBMbG9u+YZuhqKgIO44r7ul8eBedD+/ijefj1ltvJTg4mMcff5yjjjrK7jge443nIpDpfNjL4XAw4txYSg70ovLiDOg3hJ/KN1P54XNYlsXKlSt59dVX7Y7ZLP4w7MfdWP96WnLV3zncJwUY1YIMJUB0C/YXEREf8PXXX/OXv/yF22+/PaAKfxGpb/6ChTWF/4XToP9QCAqG/kPZd9ksrD5H89prr7Fxo7vGk97H54t/Z2FfgftuPZHObS2xBBjt7g2Dc3Xf5qwbICIifiA9PZ2uXbty//332x1FRGy0OGcJlcPGQsP1PSyL6gvuwbKCmDFjhj3hWsjni3+nDbif+BtFTa/+ZrEsKwtIqVv4O1f/danAffefqEbuFxERH/Xvf/+bl156ibvvvpsBAwbYHUdEbLRzxzbo18iY/kEnAoa33nqL9957z6O5WsNfiv98albcbSgOyGvOEzgn5ua5Wa237mQJd9uhpitQw5WARUTEizkcDtLSMwiPGExQcDDhEYNJS8/A4XAAkJaWRo8ePZgyZYrNSUXEbn0HRED5Zvcbd26m74AIBgwYwIwZM/D2Zjp+Ufw7J9VGOjv/ALVX7HcZY5bX3dc5dGdjg/viqLl67+rfH21ZVpzz/rqTeAsbtvW0LCsZKG1knQAREfFCrsl72fnvUj5qJubutygfNZPsZUWMODeWtWvX8vLLLzNp0iT69DlsR2kR8XMTk8YTUrwUGhb2xhBSnMfECYlMnz6dd999l8LCZg86sYU/dfsZBUx1tt0EiDLGuPs0oIJDJwkXOP+b7Gb/2qLeGFNsWZZreBDULCBW0shxRETES9WbvOcaw9t/KJUXTadkVSa33Horffr0CajWniLSuCmTJ7Fi5WuUrMqkMnos9B0COzcTUpxHVJc9TJk8ic6dO/Pggw8yffp04uLisBrOD/ASflP8O8fppzRjP3ftOJt9dpxrCWh8v4iID1ucs4TKUTPdTt6rHDycb1a9Q1ZWFj179rQnoIh4ldDQUNatLWL+goU8kZtJ+fat9Bt4JLcn3saUyZNq12BITU3ltttu47XXXuOKK66wN3Qj/GLYj4iISEs0OXnvs7cBuOOOOzyYSES8XWhoKOlpqZRt3cKa1asp27qF9LTUeouv3XjjjRx77LHMnDmTgwcP2pi2cSr+RUQk4DQ6eW/LR/D9J3TvGUb37t09H0xEfFqnTp1IT0/n008/JS+vWT1nPE7Fv4iIBBy3k/eMgX88DZ26cs9dd9oXTkR82tixY/ntb39LWloaVVVVdsc5hIp/EREJOFMmTyKqyx5CVmVC2SaoroKPX4EfPmfQwAFMvf+wU8hERNwKCgpi9uzZbNq0iWeffdbuOIdQ8S8iIgHHNXkveUws4WsysR69hE7v5dK7d28++9fH9cbwioi01P/93/9x2mmn8cgjj3hd338V/yIiEpDqTt5bvjyfqgMHeOSRR+jdu7fd0UTEx1mWxXPPPUdRUZHXtfz0m1afIiIirVFdXU1qairHH38811xzjd1xRMRP/PrXv7Y7glsq/kVEJKDl5eXx2WefsWzZMoKDg+2OIyLSoTTsR0REAlZVVRVpaWmccsop/OEPf7A7johIh9OVfxERCVjPPfcc33zzDa+++ipBQboeJiL+T7/pREQkIB04cICMjAyGDx/O7373O7vjiIifcDgcpKVnEB4xmKDgYMIjBpOWnoHD4bA7GqAr/yIiEqCeeeYZ/vOf/7B48WKv68YhIr7J4XAw4txYSg70onLUTOg3hPLyzWQvW8qKla+xbm2R7a2EdeVfREQCzr59+5gzZw4jRozg4osvtjuOiPiJ+QsW1hT+F06D/kMhKBj6D6XyoumU7O/J/AUL7Y6o4l9ERALP008/zZYtW8jIyNBVfxFpN4tzllA5bCw0/L1iWVRGj+WJ3D/ZE6wOFf8iIhJQKisryczM5JxzziEuLs7uOCLiR3bu2Ab9hrjf2HcIO7dv82wgN1T8i4iIXzncZLslS5awdetWXfUXkXbXd0AElG92v3HnZvoOjPBsIDdU/IuIiN9wTbbLzn+X8lEzMXe/RfmomWQvK2LEubH8+OOPzJ07l9jYWEaOHGl3XBHxMxOTxhNSvBSMqb/BGEKK87g98TZ7gtWh4l9ERPzG4SbbXXf99Wzfvp309HS7o4qIH5oyeRJRXfYQsioTyjZBdRWUbSJkVSZRXfYwZfIkuyOq+BcREf/R5GS7k6+koKCAuLg4zjvvPHsCiohfCw0NZd3aIpLHxBK+JpOgRZcSviaT5DGxXtHmE9TnX0RE/EiTk+2+/xRz8KCu+otIhwoNDSU9LZX0tFS7o7ilK/8iIuI3Gp1st/+/sP4lOncN4ayzzvJ8MBERL6HiX0RE/Eajk+0+egX2/cwN119vTzARES+h4l9ERPyG28l23/8L3n+W0J69ePgh+1fXFBGxk8b8i4iI33BNtpu/YCFP5Gayc/s2Qrp3578Hq3jzb3/1isl2IiJ20pV/ERHxK67JdmVbt7BzZzmdgyyuvPJKzjnnHLujiYjYTsW/iIj4rYceeoiffvpJHX5ERJxU/IuIiF/atWsXDz30EKNHj+bkk0+2O46IiFdQ8S8iIn5pwYIFOBwO0tLS7I4iIuI1VPyLiIjfKS8v55FHHmHcuHGcdNJJdscREfEaKv5FRMTvzJ8/n19++YXUVO9cYVNExC4q/kVExK+UlZXx2GOPcfXVV3P88cfbHUdExKuo+BcREb+SnZ1NZWWlrvqLiLih4l9ERPzG9u3bWbx4Mddddx3HHXec3XFERLyO36zwa1lWGDAVKHHeFWWMSWnvx7f1OCIi0nGys7PZv38/M2fOtDuKiIhX8pviH8gHkowxpQCWZUVallVgjIlv58e39TgiItIBfvjhB5544gmuv/56hg4danccERGv1OJhP5ZlXWVZ1gUdEaa1LMtKAEpdBTlAneI8ob0e39bjiIhI6zkcDtLSMwiPGExQcDDhEYNJS8/A4XAA8MADD3DgwAFmzJhhc1IREe/VmjH/E4BCy7KqLctab1nWPMuyLrAsq2d7h2uBscBGN/cXAEnt+Pi2HkdERFrB4XAw4txYsvPfpXzUTMzdb1E+aibZy4oYcW4sX3/9NTk5Odx4441ERUXZHVdExGu1uPg3xlwIRPG/Qng0UAjstixrk2VZT1iW9fv2jXlYcUCpm/tLgeHt+Pi2HkdERFph/oKFlBzoReWF06D/UAgKhv5DqbxoOiX7e3Ld9ddTXV2tq/4iIofRqm4/xpjNxpjlxpgJxpihQG9q3gysoeYK+ArnJwOLO/oTAecE3DBgl5vNFc5tbX58W48jIiKttzhnCZXDxoJl1d9gWVT+Op71H37IzTffzJAhQ+wJKCLiI9ql1acx5ifnm4EkIIaaNwBjgKHAZsuyTmmP4zSiz+F2cBbubX18W49ji507d1JYWGh3DBGRNtm5Yxv0a6Sw/+afAEyfPt2DiUREfFO7d/sxxnxkWdYYY8wSaj4BiAOWW5YVb4z5T3sfz9tZlpUIJAIMGDCAoqIijx7/6aef5vnnnycqKkpXxLyEw+Hw+PeBNE7nw7s0dj569gnnp/LNNUN+6tqzAz57ky4hR7B582Y2b97smaABQD8b3kXnw7v48vlocfFvWdYDQDJggOXAUmPMK43tb4wptCzrNCCFmv74AcUYkwvkAgwfPtzExsZ69Pi//e1vWbFiBW+88Qb5+fkePba4V1RUhKe/D6RxOh/epbHzcc+dE8letpTKi6bXH/rzwQtgDnL7hAk6j+1MPxveRefDu/jy+WjtsJ8o4HZqxvqvqNP5J8+yrFVAZN2djTEVuJ8o2x7cjcGvx3n8tj6+rcexRd++ffnDH/7A8uXL+eSTT+yOIyLSKlMmTyKqyx5CVmVC2SaoroLS9+HTN+jdtx9zZmfYHVFExCe0ZcJvrrPzT29q3ghsdv77I2PMWKj5lMCyrMXOh3VI7zVnwV1BgzccTpHObW1+fFuPY6fRo0fTs2dP0tPT7Y4iItIqoaGhrFtbRPKYWMLXZBK06FJC3p5HcHAQ77+3ltDQULsjioj4hNa0+rzfsqwplmXd5vz/n5xvBMYYYy40xtxfZ/dSYJxlWU8Aee2U2Z0NuJ+QG0VNG9L2enxbj2OLHj16cN999/HKK6/w0Ucf2R1HRKRVQkNDSU9LpWzrFjZt+pqq/fu4/fbb+fWvf213NBERn9HaK//zgdWWZZ16mP1yjTF9jDG3G2M6surMB+Ld3B9H8950NPfxbT2Obe6991569erFrFmz7I4iItJmc+bMITg4mKlTA24qmYhIm7S61adz6M/H7Zil1ZyTaiMty6odkmNZVjSwyxizvO6+lmWVWJa1sTWPb8lx7OZwOEhLzyA8YjAXXDCKY0/8LcOio3nttdfYuNHdIsUiIr7hm2++4dlnn2XChAlERETYHUdExKe0e6tPG40CplqWVeL8/yhjjLur9BW4n7zb3Mc3dz/bOBwORpwbW7Ma5qiZ0G8I5eWbeX/98wQHBzN9+nTeeustu2OKiLTKnDlz6Ny5MykpKXZHERHxOX5T/Dsn5B72L4ExJqaNj2/Wfnaav2BhTeF/4bT/tcTrP5R9l6bR6fnxrFq1ig8//JDTTz/d3qAiIi20adMmnnvuOe655x4GDRpkdxwREZ/TLiv8indZnLOEymFj6/fCBrAsqi64BysoSGP/RcQnZWRk0LVrV131FxFpJRX/fmjnjm3Qr5HVfAeeAMbw5ptv8v7773s2mIhIG3zxxRe8+OKL3HnnnQwYMMDuOCIiPknFvx/qOyACyhtZ4n7nZvoOiKBfv366+i8iPiUjI4MjjjiCKVOm2B1FRMRnqfj3QxOTxhNSvBSMqb/BGEKK85g4IZHk5GTefvtt/vGPf9gTUkSkBT777DPy8vK46667CA8PtzuOiIjPUvHvh6ZMnkRUlz2ErMqEsk1QXQVlmwhZlUlUlz1MmTyJiRMn0r9/f9LS0uyOKyJyWOnp6XTv3p0//vGPdkcREfFpKv79UGhoKOvWFpE8JpbwNZlYiy4hfE0myWNiWbe2iNDQULp3705KSgqrV6/m73//u92RRUQa9a9//Yv8/Hzuvfde+vbta3ccERGfpuLfT4WGhpKelkrZ1i2sWb2asq1bSE9LJTQ0tHafCRMmMHDgQF39FxFbNVyUMDxiMGnpGTgcDgBmzZpFz549mTRpks1JRUR8n4r/ANatWzfuv/9+ioqKeOedd+yOIyIByLUoYXb+u5SPmom55y3KR80ke1kRI86N5R//+AevvPIK9913H71797Y7roiIz1PxH+ASExMZNGgQaWlpmIYThEVEOli9RQn7D4WgYOg/lMqLplOyvye33HorYWFh3HvvvXZHFRHxCyr+A9wRRxzBtGnTWLt2LWvWrLE7jogEmKYWJaz81Qi+/uorJk+eTFhYmC35RET8jYp/4bbbbuPII48kNTVVV/9FxKOaXJTwi9UA3H333R5MJCLi31T8CyEhIUyfPp1//vOfvP3223bHEZEA0uiihNs+h/98SPcevejZs6fng4mI+CkV/wLALbfcwtFHH62r/yLiUY0uSvj+MxDcmbvuvMOWXCIi/krFvwDQtWtXZs6cyYcffsjf/vY3u+OISIBwuyjhv9+EbzcyYEB/pk+bandEERG/ouJfat14441ERkbq6r+IeIy7RQk7Fy2ie/fufFK8sd7aJCIi0nYq/qVW586dSUtL46OPPuKVV16xO46IBIi6ixIuXLCAA/v3MXv2bAYMGGB3NBERv6PiX+q55pprOO6440hLS+PgwYN2xxGRAPPMM88wcOBAJkyYYHcUERG/pOI/wDkcDtLSMwiPGExQcDCDjh7Cb357Mv/+97/Jz8+3O56IBJB33nmHTz75hKlTp3LEEUfYHUdExC91sjuA2MfhcDDi3Nia1TVHzYR+Qygv38zfNr5E15AQUlNTSUhIIDg42O6oIuLnjDGkpqbSr18/EhMT7Y4jIuK3dOU/gM1fsLCm8L9wGvQfCkHB0H8o+y6ewcF+UXz99de8+OKLdscUkQBQUFDAe++9xzXXXENISIjdcURE/JaK/wC2OGcJlcPGgmXV32BZHDj/DoI7dSY9PZ2qqip7AopIQDDGMGPGDI4++mguu+wyu+OIiPg1Ff8BbOeObdBviPuN/aI4WF1NSUkJzz77rGeDiUhA+etf/8r69etJTU2lS5cudscREfFrKv4DWN8BEVC+2f3GnZvpOzCC0047jYyMDPbv3+/ZcCISEA4ePMjMmTMZOnQoN9xwg91xRET8nor/ADYxaTwhxUuh4YJexhBSnMfEpPFkZGTw7bff8vTTT9sTUkT82ooVK/jkk0+YNWsWnTt3tjuOiIjfU/EfwKZMnkRUlz2ErMqEsk1QXQVlmwhZlUlUlz1MmTyJiy66iLPOOovMzEwqKyvtjiwifqS6upq0tDROPPFExo0bZ3ccEZGAoOI/gIWGhrJubRHJY2IJX5NJ0KJLCV+TSfKYWNatLSI0NBTLssjIyOD7779nyZIldkcWER/ScB2R8IjBpKVn4HA4AHjppZf44osvSE9PV0thEREPUZ//ABcaGkp6WirpaamN7nPBBRdw/vnnM3fuXG699Va6devmwYQi4osaW0cke9lSVqx8jbVrCpg1axannnoqV111ld1xRUQChq78y2FZlsXs2bPZvn07TzzxhN1xRMQHNLaOSOVF0ynZ35Mbb7qJkpISMjIyCArSnyIREU/Rb1xplnPPPZf4+HiysrJqP7IXEWlMU+uIVJ5yFX/72984/fTTufzyy+0JKCISoFT8S7NlZGTw448/8thjj9kdRUS8XJPriPzwJQerq5kzZw5WwzcHIiLSoVT8S7ONGDGCSy+9lPnz5/PTTz/ZHUdEvFij64gc2AcfPE/nLl2Ji4vzfDARkQCn4l9aZPbs2ezatYuHHnrI7igi4sUaXUfkk9fgl5+45pprdNVfRMQGKv6lRaKjo/nDH/7AggULKC8vtzuOiHgpt+uIbP0U/vEU3Xv04LFFj9odUUQkIKn4lxbLyMjgv//9L1lZWXZHEREv5W4dke5vpEH1AV5/9VVCQ0PtjigiEpD8ovi3LCvMsqwsy7ISnbcWVaWWZUVblpXsfI4Cy7JyLMsKc7NfifP5w5y3hMb29Wcnnngi1113HY899hjbtm2zO46IeCnXOiJlW7ewc2c5nYMtLr/8ckaOHGl3NBGRgOUXxT+QD+QYY3KNMblAjmVZBc15oGVZ0cBwY0y2MSbFGBMPlAC7LcuKbLB7JJAD7HbesoAsY0xFe70QXzFr1iyqnd06REQO56GHHqKiooKMjAy7o4iIBDSfL/4ty0oASo0xpa77XP92bjucOOcbhlrGmGygmJpCv64UIAaIB2KMMVF1j+vPHA4HaekZhEcMJig4mNPPPo9TTjmVJUuWUFoaEF8CEWmlnTt38tBDD5GQkMCwYcPsjiMiEtA62R2gHYwF3F3lLwCSgOWHeXySZVnFxpjCBvcXAskNdzbGFLcqpQ9zOByMODe2ZrXOUTOh3xDKyzez94O/UF1dzYwZM3jxxRftjikiXmr+/Pk4HA5mzZpldxQRkYDn81f+gTjA3aXnUmB4M5+j4fAeqWP+goU1hf+F06D/UAgKhv5D2Xd5BkF9jmLp0qV8/vnndscUES+0fft2Fi1axDXXXMNJJ51kdxwRkYDn08W/c6JtGLDLzeYK57YmOYfu5LrZFEnN0B93x422LCtgVqdZnLOEymFjoWFPbsuieuQ9AKSmptqQTES8XWZmJvv379dVfxERL+HTxT/Q53A7tKYTj/MxCcC8BptOc84jqDDGFLo6DLX0+X3Nzh3boN8Q9xuP/C3GwIoVK9i4caNng4mIV9u8eTM5OTnceuutDB061O44IiICWKbh6os+xNmNp4SaybfFDbbFUTPuv3dLu/FYlpUD9DHGjG74nA3nBliWVQIkuZkz4NqeCCQCDBgwIGbp0qUtidIuHA5Hm3pqXzl6HD9dnFEz5Kehsk30ejMVc6CS448/Xr3/m6Gt50Pal85Hx5k3bx5FRUU8//zzhIeHN+sxOh/eQ+fCu+h8eBdvPx8jR47caIxxO/zdHyb8tivnm4bhxpiYhtsaKfCXU9Py85D9nY/JBXIBhg8fbmJjY9svbDMVFRXRluPec+dEspctpfKi6fWH/hhDSHEe99x1B92PCCElJYVOnTpxzjnntD20H2vr+ZD2pfPRMT777DMKCgr44x//yOjRow//ACedD++hc+FddD68iy+fD9uH/TgXzmrJre7YEndj/etpyVV/53CfFGBUC15CCRDdgv19zpTJk4jqsoeQVZlQtgmqq6BsEyGrMonqsocpkydx5513MnDgQKZNm4Yvf5okIu1jxowZ9OjRg5SUFLujiIhIHbYX/84Jty25xdR5bAU1E3vddeuJdG5riSXAaHdvGJxvPJqzboDfCQ0NZd3aIpLHxBK+JpOgRZcSviaT5DGxrFtbRGhoKN26dWPGjBmsXbuWt99+2+7IItKBGq77ER4xmLT0DBwOBwAffPABK1euZMqUKfTt29fmtCIiUpftxX872ID7ib9R1PTqbxbLsrKAlLqFv3P1X5cK3Hf/iWrkfr8SGhpKeloqZVu3UF1dRdnWLaSnpdYb7zZ+/Hh+9atfMX36dF39F/FTrnU/svPfpXzUTMzdb1E+aibZy4oYcW4sDoeDadOm0b9/f+69916744qISAP+UPznU7PibkNxQF5znsA5KTfPzWq9dSdKuNsONV2BGq4EHJC6dOnCrFmz2LhxI6+88ordcUSkAzS27kflRdMp2d+T2yfewZo1a5g+fbpXT4YTEQlUPl/8OyfURjo7/wC1V+x3GWPqre7rZs6Aa4JvlOtxrh7+zvvrTuItbNjW07KsZKC0kXUCAtJ1113Hr3/9a2bOnEl1dbXdcUSknTW17kflsDEsXbqUX/3qVyQlJdkTUEREmuQv3X5GAVOdbTcBoowx7j4NqODQScIFzv8mu9m/tqg3xhRbluUaHgQ1C4iVNHKcgNWpUydmz57NmDFjeOmll7juuuvsjiQi7ajJdT8qtlF1oGZBr65du3o2mIiINItfFP/OcfqHbSnRSPtOy92+jTy+mAAY399Wf/jDHzj11FNJS0tjzJgxdOnSxe5IItJO+g6IoLx886HrfhyshveWENypE9dff7094URE5LB8ftiPeJ+goCDmzp1LaWkpS5YssTuOiLSjiUnjCSleCg0n9X9eAD9t56qr/kBwcLA94URE5LBU/EuHuPjiizn//PPJyPhf+z8R8X1u1/344XN451GO6NaNp/6kN/wiIt5Mxb90CMuyyMrKoqysjIULF9odR0Taibt1P7r/bSYc2Efe0qX06NHD7ogiItIEFf/SYc444wyuuuoq5s+fT1lZmd1xRKSd1F3346efKujWtTMXXHABv/vd7+yOJiIih6HiXzrU3Llz+eWXX5gzZ47dUUSkAzz88MP8+OOPzJ071+4oIiLSDCr+pV05HA7S0jMIjxhMUHAw54yM45RTT+XJJ5+ktNTdGmki4qt27tzJ/PnzueKKKzjjjDPsjiMiIs3gF60+xTs4HA5GnBtbs/rnqJnQbwjl5ZvZ+8FfqKqqYurUqeTlNWvRZRHxAVlZWezdu1ef7ImI+BBd+Zd2M3/BwprC/8JpNT3Ag4Kh/1D2XZ5BUJ/BLFu2jI8++sjumCLSDr7//nsWLVrEddddx29+8xu744iISDOp+Jd2szhnCZXDxoLVYN00y6J61L1YVhBTp061J5yItKu0tDQOHjxIRkaG3VFERKQFVPxLu9m5Yxv0G+J+46CTAMOqVatYs2aNR3OJSPv69NNPeeaZZ7jrrrs45phj7I4jIiItoOJf2k3fARFQvtn9xp2b6TvwSAYPHkxKSgqm4eqgIuIz7r//fnr27Mm0adPsjiIiIi2k4l/azcSk8YQUL4WGhb0xhBTnMTFpPBkZGWzYsIHly5fbE1JE2mTNmjW88cYbTJ8+nT59+tgdR0REWkjFv7SbKZMnEdVlDyGrMqFsE1RXQdkmQlZlEtVlD1MmT+L666/npJNOYvr06Rw4cMDuyCLSAgcPHiQ5OZmjjz6aO++80+44IiLSCir+pd2Ehoaybm0RyWNiCV+TSdCiSwlfk0nymFjWrS0iNDSU4OBg5s2bx6ZNm3jqqafsjiwidTRcpyM8YjBp6Rk4HA4A8vLy2LhxI3PmzCEkJMTmtCIi0hrq8y/tKjQ0lPS0VNLTUhvd5/LLL+ecc84hPT2d66+/nu7du3swoYi409g6HdnLlrJi5Wu8W7iK6dOnc8opp3DttdfaHVdERFpJV/7F4yzLIisri+3bt/Pwww/bHUdEaHydjsqLplOyvyfXXHsdmzdvZv78+QQF6U+HiIiv0m9wscVZZ53FFVdcQVZWFuXl5XbHEQl4Ta3TUfmb31Hw9ttceOGFxMfH2xNQRETahYp/sc3cuXP5+eefmTt3rt1RRAJek+t0lKzDmINkZWV5NpSIiLQ7Ff9imxNPPJGbbrqJxx9/nNLSUrvjiAS0Rtfp2LMDPnqZrkd049RTT/V4LhERaV8q/sVWGRkZdOrUifvvv9/uKCIBrdF1Ov75DJiDJCUl2ZJLRETal4p/sdWRRx5JcnIy+fn5/OMf/7A7jkjAcrtOx1fvwOdv0zc8nMzZGXZHFBGRdqDiX2z3xz/+kYiICCZNmsTBgwftjiMSkNyt09G58EFCQo7gow3rCQ0NtTuiiIi0AxX/Yrvu3buTmZnJhx9+SF5ent1xRAKWa52Osq1bWLXqLQ7sqyQzcw6DBw+2O5qIiLQTFf9ii4Yrif5x6gwGDhxISkoKv/zyi93xRALawYMHSU5O5le/+hV33HGH3XFERKQdaYVf8Th3K4nuLN9M57VPcODbj8nOziYtLc3umCIB66WXXuKjjz7i+eefp2vXrnbHERGRdqQr/+Jxja0keuCqBwkK7cvcufMoKyuzO6ZIQKqsrGT69OkMGzaMq6++2u44IiLSzlT8i8c1tZLowZF3sX//Pl35F7HJQw89xLfffsuDDz5IUJD+RIiI+Bv9ZhePa3Il0cgzAYvc3Fw+++wzj+YSCXQ//PADmZmZXHnllVxwwQV2xxERkQ6g4l88rtGVRAF2bqbvgEH06NGDKVOmeDaYSICbPn06+/fvZ/78+XZHERGRDqLiXzyu0ZVEjSGkOI87bk9i5syZvPnmm7z99tv2hBQJMBs3buSZZ57h3nvvZejQoXbHERGRDqLiXzzO7UqiZZsIWZVJVJc9TJk8iTvvvJPIyEgmT55MdXW13ZFF/JoxhnvuuYfw8HBmzJhhdxwREelAKv7F49ytJBq+JpPkMbGsW1tEaGgoXbt2JSsri3//+988/fTTdkcW8WkN19UIjxhMWnoGDocDgPz8fP7xj38wZ84cevbsaXNaERHpSOrzL7ZwrSSanpba6D5/+MMfOPvss5kxYwbjxo2jR48eHkwo4h/cratRXr6Z7GVLWbHyNd55+02mTJnCKaecwi233GJ3XBER6WB+ceXfsqwwy7KyLMtKdN6yWvj4Eufjwpy3BMuycizLCmvP40jLWJbFwoULKSsrIytLX2qR1mhsXY3Ki6ZTsr8nY8eNY8uWLTz88MMEBwfbHVdERDqYXxT/QD6QY4zJNcbkAjmWZRW04PGRQA6w23nLArKMMRXtfBxpodNPP51rrrmGBQsW8N1339kdR8TnNLWuRuXxF/HOO+9w1VVXERsba0s+ERHxLJ8v/i3LSgBKjTGlrvtc/3Zua44UIAaIB2KMMVF1n68djyOtMHfuXIwxTJ061e4oIj6nyXU1vlgNxqi1p4hIAPH54h8YC2x0c38BkNTcJzHGFBtjCo0xxR15HGm5X/3qV0yaNIkXXniBf/7zn3bHEfEpja6rsf1L+KKAI0J7EBkZ6flgIiJiC38o/uOAUjf3lwLDffA44sa0adM48sgjufPOO9X6U6QF3K6rYQy8sxiCu3D3XXfZF05ERDzOp4t/54TcMGCXm80Vzm0teb5oy7LiOvo40nKhoaE8+OCDfPTRR/zpT3+yO46Iz3C7rsb6l+CHz4iIGMSMaRpOJyISSHy6+Af6HG6Hhh17GnGac9x+hTGm0NXRpwOOI20wduxYzj//fKZNm8bOnTvtjiPiExquq2E9eglB7z/DwIGD+PzTTwgNDbU7ooiIeJBl6n4U7GMsy4oESqiZpFvcYFscNePxe7vp2tPweeKMMYUN7isBkpxvBlp9HOebiESAAQMGxCxdurRlL7IdOBwOv/kDX1payvjx47n88su577777I7TKv50PvxBoJ2PZ599lj//+c888sgjnHzyyXbHOUSgnQ9vpnPhXXQ+vIu3n4+RI0duNMa4HZau4r/x584C4owxMe11nOHDh5sNGza0NEqbFRUV+WwbP4fDwfwFC1mcs4SdO7bRd0AEkUcfyfoPP2TDhg1ER0fbHbHFfPl8+KNAOh9bt27luOOO47LLLmPZsmV2x3ErkM6Ht9O58C46H97F28+HZVmNFv+2r/DrvMLeEhXGmBjnv92Nwa+nNYW/UwmQ7IHjSCMaW5l07/rnCAoOZuLEibz//vtYDfuXi4hb999/P9XV1WRnZ9sdRUREbGL7mH9nT/2W3GLqPLaCmgm37vrURTq3Ncm5um+Tffrb4zjSco2tTLrv0llY4UP54IMPeP755+2OKeIT/vGPf/D8888zefJkjjnmGLvjiIiITWwv/tvBBtxPyI0CCt3c31AF4K63f1SD+9t6HGmhplYmrbrgbjp17sKUKVPYs2ePPQFFfERVVRW33347Rx99NNOmTbM7joiI2Mgfiv98albmbSgOyGvG4/MarubrlADktONxpIWaXJm0XxTVVVWUlZWRnp7u2WAiPmbRokV8+umnPPLII3Tv3t3uOCIiYiOfL/6NMblApHNSLlDTrx/YZYxZXndf5xCfhqv0FjZo64llWclAqfO5W3wcaR+NrkwKsHMz/QYdya233sqjjz7K559/7tlwIj5i69atpKamctlll3HFFVfYHUdERGzm88W/0yggybKsRGchP9YY4+4qfQUNJu86u/dscPb2z7IsK8d5v7vHN/c40g7crkwKYAwhxXncnngbc+fOJTQ0lLvvvhtf7lwl0lEmTZpEVVUVjz76qCbHi4iIfxT/xpgKY0yKMSbXeUtpZL8Yd8W6MabY+fgUY0ySMcZtK4zmHkfah9uVScs2EbIqk6gue5gyeRLh4eHMnj2b1atXs2LFCrsji3iUw+EgLT2D8IjBBAUHEx4xmLT0DBwOBwAFBQUsW7aMadOmERnprl+BiIgEGr8o/sU/NVyZNGjRpYSvySR5TCzr1hbVLq4xYcIETj75ZCZNmsR///tfWzOLeIqrFW52/ruUj5qJufstykfNJHtZESPOjWXnzp3ccccdHHvssUyZMsXuuCIi4iVU/ItXCw0NJT0tlbKtW6iurqJs6xbS01LrrarXqVMnFi1axHfffce8efNsTCviOY21wq28aDol+3uSMHoMmzZt4rHHHiMkJMTuuCIi4iVU/ItfOO+887jmmmuYP38+JSUtXTdOxPc01Qq38rhRFL3zDqNHj+bCCy+0J6CIiHglFf/iN+bPn0/nzp259957NflX/F6jrXCNgY9fAwwPPfSQx3OJiIh3U/EvfiMiIoK0tDT++te/8sorr9gdR6RDNdoKt+Sf8J8P6d4zjCOPPNLzwURExKup+Be/cu+993Lqqady5513UlFRYXcckQ7jthXugV/gncewunbn3nvuti+ciIh4LRX/4lc6derEkiVL2LFjB1OnTrU7jkiHcdsKd80i2FvGrwYfxf3J6vAjIiKHUvEvfmf48OHcc889PPnkk7z33nt2xxHpEA1b4VqPXgKfreKUU0/l04821OuIJSIi4qLiX/xCw8WOnn0xj169enHbbbexb98+u+OJdAhXK9wd339LbOz5hIWF8faqVSr8RUSkUZ3sDiDSVq7FjkoO9KJy1EzoN4Sd5Zvp/O7jfPXVv5g9ezZz5syxO6ZIh3nppZd45513eOKJJ+jfv7/dcURExIvpyr/4vMYWOzqQsJCgnv154IEH+PLLL+2OKdIhdu/ezaRJkzjttNMYP3683XFERMTLqfgXn9fUYkcH4//IwYMHSUxM5ODBg/YEFOlAkyZNory8nCeffJLg4GC744iIiJdT8S8+r9HFjgCOOgWAtWvX8tRTT3kwlUjHe+utt3jmmWdISUkhOjra7jgiIuIDVPyLz2t0sSOAnZvpO/BIYmNjmTJlCj/88INnw4l0kD179pCYmMgJJ5zAzJkz7Y4jIiI+QsW/+Dy3ix0BGENIcR4Tk8aTk5NDZWUl99xzjz0hRVqgYfeq8IjBpKVn4HA4avdJSUnh+++/5+mnnyYkJMTGtCIi4ktU/IvPc7vYUdkmQlZlEtVlD1MmT+K4445jxowZ5Ofn8/rrr9sdWaRRru5V2fnvUj5qJubutygfNZPsZUWMODcWh8NBUVERTz75JPfddx8jRoywO7KIiPgQFf/i8xoudhS06FLC12SSPCaWdWuLanueJycnc9JJJ3HHHXewd+9eWzOLNKax7lWVF02nZH9PMuc9wK233kpUVBSzZ8+2O66IiPgYFf/iF1yLHZVt3UJ1dRVlW7eQnpZab7GjLl26sGTJEr7//ntmzJhhY1qRxjXVvaoyeiyPPPIIpaWlPPXUU3Tr1s2ekCIi4rNU/EtAOfPMM7n99ttZtGgRH374od1xRA7RZPeq/b/wy88Obr/9ds4//3zPBhMREb+g4l8Czrx584iIiOC2225j3759dscRqafR7lVV++GtLIKCg8nKyvJ8MBER8Qsq/iXg9OzZkyeffJJPP/2UtLQ0u+OI1NNo96r3n4U92xl39TX06NHDnnAiIuLzVPxLQLr88su57bbbyM7O5r333rM7jkgtt92rviiE9UsJ69OXnCcW2x1RRER8mIp/CVgLFy7kmGOO4cYbb1T3H/EaDbtXWY9eQnDBg4SGhvKvjz+qN4ldRESkpVT8S0Cpu3hSr7Awdu1xUFq6mbvvvtvuaCK16navmjUrjeqqA7zwwvMMHjzY7mgiIuLjOtkdQMRTXIsnlRzoReWomdBvCD+Vbyb4zUyeeeYZLrvsMhISEuyOKVLr3//+N3PmzGHcuHH83//9n91xRETED+jKvwSMxhZPqr7mSayu3bnpppspLy+3O6YIAAcOHODmm2+mV69ePProo3bHERERP6HiXwJGo4snde6KibuPn392MGHCBEzDLisiNkhPT2fDhg0sXryY8PBwu+OIiIifUPEvAaPJxZOGnotlWaxYsYIXXnjBs8FEGnj33XeZO3cut9xyC6NHj7Y7joiI+BEV/xIwGl08CWDnZvoOPJKzzz6bO++8k++++86z4UScdu3axXXXXcfQoUN55JFH7I4jIiJ+RsW/BIxGF08yhpDiPCYmjecvf/kLVVVV3HzzzRw8eNCeoOK36nabCgoOJjxiMGnpGTgcDgCMMYwfP54dO3bw0ksvqa2niIi0OxX/EjDcLp5UtomQVZlEddlTsz0qioceeojVq1fz2GOP2R1Z/Iir21R2/ruUj5qJufstykfNJHtZESPOjcXhcPDUU0/x8ssvM2fOHGJiYuyOLCIifkjFvwSMhosnBS26lPA1mSSPiWXd2qLaq6y33XYbl112GSkpKXz55Ze2Zhb/0Vi3qcqLplOyvycp90/lnnvu4YILLuCPf/yj3XFFRMRPqfiXgFJ38aTq6irKtm4hPS213vAKy7L405/+RPfu3bn++us5cOCAjYnFXzTabcqyqDzlKnJzczjiiCN49tlnCQrSr2YREekYAf8XxrKsSLsziPcZOHAgOTk5bNiwgczMTLvjiB9ostvUV0VUHTjAU089xZFHHunZYCIiElD8YoVfy7LCgKlAifOuKGNMSjMfnmRZVjSQA1QAuxruYIwpdh6nBMgCljk3xQHxQIoxpqKV8cVL/eEPf+CGG25g9uzZnH/++YwcOdLuSOLD+g6IoLx8c82Qn7r+swGKVxDSrTtXXHGFPeFERCRg+MuV/3wgxxiTa4zJBXIsyypo5mMjqSni84ECYGODW1KDfXOA3c5bFpClwt9/Pf744/z6179m3LhxbNu2ze444sPcdpv672546wGsrt247957bcsmIiKBw+eLf8uyEoBSY0yp6z7Xv53bDqfUGGO5u1FzRb9u8Z8CxFBztT/GGBNV97jif0JDQ1m+fDkOh4Nx48ZRVVVldyTxUYd0m6o6AK+nwy8/EXnMr5g29X67I4qISADw+eIfGEvNFfqGCqh/1b4xJe7utCwrDihueL8xptgYU+gaCiT+qW4/9t/89rcEdQlh7dq1TJkyxe5o4qMadpuyFl0CWz/loosu4uMNH6qnv4iIeIQ/FP9xgLur76XA8MM92DlMyJ1oY0xhW4KJb3LXj91x+TyCe0fw8MMPk5eXZ3dE8VGublNrVr1Bly5duOSSS3jzjb+p8BcREY/x6eLfOdE3DDeTdKmZvBvWyudNNsZkN7E92vnJgPihxvqxV1/3J6yQHtx0002Ulmq0l7SOawhZWFgYzzzzDFbD1p8iIiIdyKeLf6DP4XZwvkFoNmfnn8Yqu9Oc8wgqjDGFlmVlWZaV2JLnF+/XaD/2zl0xF6ewb98+Ro8eTWVlpT0BxWcZY7j55pv54osveO655+jfv7/dkUREJMBYpm7nCR/j7NFfQs3k2+IG2+KoGfffuyXdeCzLyjfGjG5kW1zDoUDO9p9JjQ0Rcr45SAQYMGBAzNKlS5sbpd04HA4NK2iBCy4YhbnnrZor/g1VV8GjF4Mx/O53v2PSpEktfn6dD+/iyfPxwgsv8Kc//YmkpCTGjRvnkWP6Gv18eA+dC++i8+FdvP18jBw5cqMxxu3wd7/o899eDrfgVyMF/nJqWn7GNPKYXCAXYPjw4SY2NraNKVuuqKgIO47rq/oObKQfO8DOzYQPOopbb7iWBx54gNGjR3P99de36Pl1PrxLe5wPh8PB/AULWZyzhJ07ttF3QAQTk8YzZfKk2j8Ob775Jk899RTjxo3jiSee0HCfRujnw3voXHgXnQ/v4svnw/bi33nlvCUqjDGuQtvdWP96WtiDPwlY38I8JUByCx8jXmxi0niyly2l8qLp9Yf+GENIcR63J97GzOnTeP/995kwYQLDhg3jN7/5jX2BxVauCeIlB3pROWom9BtCeflmspctZcXK11i3togffviBq6++mlNOOYWnnnpKhb+IiNjG9jH/zl75LbnF1HlsBTUTe91dsY90bmuJBBoZ729ZVkkz1w0QH3dIP/bqKijbRMiqTKK67GHK5El06tSJpUuX0rNnTxISEti7d6/dscUmjU0Qr7xoOiX7ezJn7jyuvPJKOnXqxCuvvEK3bt3sjiwiIgHM9uK/HWzA/cTfKKDZrTqdE4ObesNQgZu+/87jqOe/H2nYjz1o0aWEr8kkeUws69YW1Q7jGDhwIEuXLmXTpk3cdttt+PL8GWm9RieIWxaVw0bz0EMP8dVXX7Fs2TKOOeYYWzKKiIi4+EPxn0/NirsNxQEtacju+vSgsaFEeY2s5psA5LTgOOIDXP3Yy7Zuobq6irKtW0hPSz1kcs/5559PZmYmy5Yt47HHHrMprdhp545t0G+I+40l69hf+QsPPvggF1xwgWeDiYiIuOHzxb9zQm1k3cm6znadu4wxy+vu6xy64241YPhf8V/RyPbChm09LctKBkqbWChMAkBycjK/+93vuO+++3j77bftjiMe1ndABJRvPnRDyT9h3bN0PaIb99xzj+eDiYiIuOHzxb/TKCDJsqxEZ4E+1hjj7tOAChq/sl/qvLnd7mwlusHZ2z/Lsqwc5/3ujiMBJCgoiBdeeIHf/OY3JCQk8Mknn9gdSTxoYtJ4QoqXQt1hX7u2wJvzsEJ6MGnSJE3wFRERr2F7t5/24Jz4m9KM/dy243RuK6Zm/H5Tjy9G4/uljrotHsu3byXICuLcc89jw4b1HHfccXbHEw+YMnkSK1a+RsmqTCqjx0Jof3j5fqg+wLFDo5h2/2F/NYmIiHiMv1z5F/E4V4vH7Px3KR81E+5ZxcFLprL355855ZRT+eGHH+yOKB5Qd4J4v9Vz4MmrYM8Obrz+OjZ+8E+vXgRGREQCj4p/kVZy2+Lx1yPhijlUVlZy7nnnUVVVZXdM8QDXBPHbb7sZgMcee4xn/vy0Cn8REfE6Kv5FWqnRFo9DToczrqHkm2+444471AI0QOTm5jJ79mxuueUWJk6caHccERERt1T8i7RSky0eR9wAWOTm5pKdne3RXOJ5y5cvZ8KECVx22WU8+eSTmuArIiJeS8W/SCs12uIRYOdm+g06knHjxnH//feTl9eSJSfEWzgcDtLSMwiPGExQcDDhEYNJS8/A4XDU7lNQUMA111zDWWedxbJly+jcubONiUVERJqm4l+kldy2eAQwhpDiPCYmjefPf/4z55xzDjfccAPvvfeePUGlVRpO6DZ3v0X5qJlkLytixLmxOBwOPvzwQ37/+99zwgkn8Ne//pVu3brZHVtERKRJKv5FWmnK5ElEddlDyKpMKNsE1VVQtomQVZlEddnDlMmTCAkJYeXKlRxzzDFcccUVfP3113bHlmZyO6G7/1AqL5pOyf6eJKfczyWXXMKAAQN46623CAsLszuyiIjIYan4F2mlui0ew9dkErToUsLXZJI8JpZ1a4tqO7307duXN954g6CgIC699FIqKirsjC3N1OiEbsui8rg4cp58ki5dulBQUMCgQYPsCSkiItJCKv5F2sDV4rFs6xaqq6so27qF9LTUQ1o8RkVF8dprr7F161amTp2qNwA+oNEJ3f/dDUVPcPBgNatWrSIyMtLz4URERFpJxb+Ih5x55pnk5eXxzTffcOGFF7J79267I0kT3E7o3vczvDwV9u4grG9/Tj75ZHvCiYiItJKKfxEPcHWNuXXCHVRVVbN+/XpOOukktmzZYnc0acQhE7qr9sOrM6G8lM4RJ3D3XXfYG1BERKQVVPyLdLCGXWO4dxWMvJMftu/g18cfz7fffmt3RHGj3oTu7V/CXzPg+0/oPOjXHNe3M1MmT7I7ooiISIup+BfpYG67xgz7PVyZSWXlPs444wx+/PFHu2NKA64J3X9MOI+ur0yB0vcJ7RnG1PFj603oFhER8SUq/kU6WKNdY4acDhfcyY4dZYwcOZIdO3bYE1Aa1alTJz7/96fs++W/pKens/en3W4ndIuIiPgKFf8iHazRrjEAv7kMy7IoLS1l5MiRbN++3bPhAlTdlXsvuGCU25V7f/rpJy655BJefvllHn74YVJTU21MLCIi0j5U/It0MLddY1x2bqbfoCN588032bJlC7GxsWzbts2zAQPMISv33nPoyr07duwgNjaW9957j+eff5577rnH7tgiIiLtQsW/SAc7pGuMizGEFOdxe+JtnH/++bz55pts3bqV2NhYtm7dak/YAHC4lXunTZ/B2Wefzddff83rr7/Otddea3dkERGRdqPiX6SD1esaU7YJqqugbBMhqzKJ6rKntmvMueeey1tvvcX27ds5//zz+e6772xO7p+aXLn3mLN47LFF7N69mzVr1nDxxRfbE1JERKSDqPgX6WCurjHJY2IJX5OJtegSwtdkkjwm9pCuMWeffTZvv/02P/74I2effTbFxcW25fZXjc7B+O4TKHwIc/Ag7733HmeccYbnw4mIiHQwFf8iHhAaGkp6WiplW7ewZvVqyrZucds1xuFw8OaqtwnqEsJ3331HTEwMf0gYXW8iqrSN2zkY37wHL6dASE/6DBjECSecYE84ERGRDqbiX8RL1J2IWnHxbEjMg/5DeXnFco6JGspPP/1kd0S/cMgcjE/fgNfTITyKrgOHcuftE+wNKCIi0oFU/It4iUMmoob2g6sfg5N/x86yHcTEDGf37t12x/R5rjkYXd+aA6sfgYIFMPB4uob1Z2i3X7Ryr4iI+DUV/yJewu1E1ODOEHcvnHEtJSXfcPrpp/P555/bltHb1e3fHxQc7LZ/f2hoKK+/nM9RB76HT14DoB97SBk3Siv3ioiI3+tkdwARqdHkYmAjbsD68CX27t3LGWecwfPPP88VV1zh2YBezjVsquRALypHzYR+Qygv30z2sqWsWPlabWG/evVqrrvuOnbv3s3jjz/OCSecwMiRI+2OLyIi4hG68i/iJZqzGNiGDRs44YQTuPLKK8nIyODgwYOeDenFDte/Pyt7PtOnTyc+Pp6wsDA++OADJk6ciNWw5aeIiIgfU/Ev4iWasxjYUUcdxd///nduuOEG0tLSSEhIYO/evfYE9jJN9u8/Lo4HHniAuXPncvPNN7NhwwZOOeUUe4KKiIjYSMW/iJdo7mJgISEhPPbYY1x40cW88sor9OzZk159ww8Z2x5oGh02tWktvDmPqgP7efHFF3nqqafo3r275wOKiIh4ARX/Il6i4WJgQYsudbsYmMPh4MzzRvL37/dD/GQIi2DPrnLmPLSY4WeeHbBvAA4ZNnVgX003n9dnQfe+9Ok/iKuvvtq2fCIiIt5Axb+IF6m7GFh1dZXbxcDqjW3/7aVw/Z/gjOs4uHcnX33xOddcex2m4dAhH9ecLj71hk3t/BZeurOmm0/MaLoOOIY7J6p/v4iIiIp/ER9zyNj2zl3h7Jvh+lzoewyvv/YqsbGxfPHFF/YGbSd1Fz8rHzUTc/dblI+aSfayIkacG1v7BmDK5En8KmgnwX+5CZ4bDz/vhJF3EvJLOUO7OtS/X0REBBX/Ij6n0bHtfX8F4xZhWRaffvopp5xyCjNnzuSXX37xfMh2dLguPvMXLOS///0vjz76KNu+3czB3dsICQnBqtxL+FevHjJsSkREJJCp+BfxMU22BN31Lf0GHcWXX37JuHHjmDNnDkcddRS9+oY3OlzG2zXZxefU0Sx8+GGOPfZYpk+fzgUXXMBnn/2bX352cPBgtdthUyIiIoFMxb+Ij2lOS9D+/fuzePFifhV1LLsdv7BnVzlmyJmU//ZqsvLeqTdcxtu5/aTDGChdB6uycFTs5uijj2bt2rWsXLmSE044wZ6gIiIiPsCvin/LshIsy8pvxePCLMvKsiwr0XnLast+Ih2puS1B5y9YyI6QwZikFXDmjfD9J/D2fPbt3MbXZf9lXla2za+keRN5D/mkY/uXkD8ZVk6H/f+lZ+++/POf/+Scc86x4RWIiIj4lk52B2gPlmXlOP8ZCfRpxVPkA0nGmFLn80VallVgjIlv5X4iHcbVEnT+goU8kZvJzu3b6DswgtsTb2PK5Em1Q1wW5yyhctTMmgnBZ94Aw0fD54Xw0cscKNvCvLmZdOkUzIQJExgwYIDHX4drIm/JgV41OfsNobx8M9nLlrJi5Wu14/QnJo0na+kL7Bs6Er4ogJJ/whFhMPIuum7/F/eOvUCr9IqIiDSTX1z5N8YkGWOSqCnOW8SyrASg1FXQO5+vtM62Fu0n4gnNaQl6yHCZzkfAKb+DG5+GK+diDh5k1qxZHH300VxzzTUkJk1o8gp8SzTniv7hJvJmzX+Q1atXU/rNJg58/R78NR22fgqnXwO/SyNkx7/VxUdERKSF/KL4b6OxwEY39xcASa3YT8QrNDox2LIgtA/hEYP56quvuPnmm8nLy2NJbg7lVSGYUfdSfnqS27kBzSnqm9ua0+1EXmNgx9dUHrSYO2c2cXFxvPrqq4wbN45rr7+Bfr17EbQhj/D3FqqLj4iISCv4xbCfNooDctzcXwoMb8V+Il5hYtJ4spctpfKi6YcU2K6JwccddxwDBkXQ+biz2TfgJPj4VShYCMC+bn34olMwo8eMYd7cuQwZMoSzY0cddphOvSv6ruO6ruivymT+goWkp6X+75OJA7/Azi01E3i/XA0VWyG4MwcPHiQ/P5/LLruMI444wvNfQBERET8U0MW/ZVlhQBiwy83mCue2Zu8n4k2mTJ7EipWvUbIqk8rosdB3COzcTEhxXr2JwYtzlrBv1MyaoTfRCbBrC2z9F3z/KQe/+4i33nyTt958k65du3KgcygHT7mipmCv3Au9j6LywmmUvD23tqivnWvQ8Ir+z7uojDiVhQ89TMWunXTq1JkDuWPhv7udO1kw+FQ4bRyERRD+j4dJSNCIOhERkfYU0MU/zZgc7Cz8m7WfMaaiHTKJtIvmTgyuNzfAsmoWC+v7Kzj5d1BdhfXoJTz33LPclpjEweAu8I+nDzlWZXAXZme8y7PP/JnyH76v+fSgc0jN7Zc9sHsL7PsZAAfw1FNP0a9vX3b8t5qDZ90MfY6GQSdAj/CaTyZWZXJ74m0e+kqJiIgEjkAv/jucZVmJQCLAgAEDKCoq8ngGh8Nhy3HFPU+fj5Hnn8fI88+rd9+GDRtq/92zTzg/lW+uufLf0M7N9OrXnyOPPJJ9v/wCia/APkfNxNu9P8KBSuftF0zxCo477ji+3/YDVcGda672/3c3dOkOx8fVFPhBQfT86EVWLs+jsrKSifdMYtuuUvZHngHdekPZJrpseImB1TsYcfppHvk66efDu+h8eA+dC++i8+FdfPl8qPjvYMaYXCAXYPjw4SY2NtbjGYqKirDjuOKet52Pe+6c2OTcgLvvuJ3Y2Fj6Doyg3PUmYWiDnvplmwjf/iGrVq0iLT2D7GVF7p9vVSb33n0nI0eOBODTc8897CcTHc3bzkeg0/nwHjoX3kXnw7v48vmwvduPZVklLby567jTWu7G8NfjHMrT3P1EfE5zFw1rzsrCLXk+aF7LUhEREWk/thf/xpioFt5i2vHYFdRM2I10sznSua3Z+4n4ItfcgOQxsYSvySRo0aWEr8k8pJVmc4v65j6fiIiIeJ6G/cAG3E/ojQIKW7GfiM9xXYFPT0ttcp/mTCBu7vOJiIiI59l+5d8L5APxbu6PA/JasZ+I39IwHREREd/mb8V/GE303Hc3Z8A5ITfSsqzIOvtFA7uMMctbup+IiIiIiLfyi2E/lmVlUVP0jwHCLMvKp2aSbo4xprjOrhW4n7w7CphqWVaJ8/+jjDHurvI3dz8REREREa/jF8W/MSbF+c+kw+zndrKwc0JvirttrdlPRERERMQb+duwHxERERERaYSKfxERERGRAKHiX0REREQkQKj4FxEREREJECr+RUREREQChIp/EREREZEAoeJfRERERCRAqPgXEREREQkQKv5FRERERAKEin8RERERkQCh4l9EREREJEBYxhi7MwQMy7J+BL614dD9gHIbjivu6Xx4F50P76Lz4T10LryLzod38fbz8StjTLi7DSr+A4BlWRuMMcPtziE1dD68i86Hd9H58B46F95F58O7+PL50LAfEREREZEAoeJfRERERCRAqPgPDLl2B5B6dD68i86Hd9H58B46F95F58O7+Oz50Jh/EREREZEAoSv/IiIiIiIBopPdAaRjWJYVBkwFSpx3RRljUuxLJJZlJQBjjTGj7c4S6CzLigbigL5ANFAKpBhjKuzMFagsy4oEkurcFUnN+Si1KZI4WZaVCJQaYwrtzhKILMsqAbKAZc674oB49PvKNnXqq53U/A0BmOdL50PDfvyUZVkFQJLrj6fzj2uOMSbe3mSBx7KsHOc/I4E+xpgYO/MEOmfhP9wYk1vnvmRq/sBGqeD0LOfvpgRjTHad+xKAJUCMzod9nEXObiBexb89LMtqWKSVUnM+9HNhA1ctRf36KgE4zZcusGrYjx9yfiOW1v3l0OCbVDzIGJNkjEkC8u3OIgDE1S38AZyFZzE1v9TFs5Ia3mGMWQ6EudsmHjUGqLA7RIBLAWKoudofY4zRBQp75XPop5Jj7QrTWir+/dNYYKOb+wvQH1ORJMuy4tzcX0jNR+riee7+eFZQ8wZAbOD8GdHVfi9gjCk2xhQaY4rtzhLIXBdPG54HY8xoX7rqDyr+/VUcNR8NNlQK+ORqdCLtLNLuAFLDGJPScCicc7hJGPq0zE6RusIsUk8SfvKGWBN+/UydP5q73GyuQFfSJMAZY6Ia2RRJzdAfsd8SIFfjzO1hWVZiw6FxYi/nXKU++pmw1XAgx/mpWCQ1ddZp+NhkX1Dx74/6HG4Hy7LCfO0bVaQjOd80JwDqxGSTOh2YTgPynOP+xcOcExp1xd97nOYcblJsjCm2LCsLKNGbM1uE8b/fT4UAlmUVAhsty4rxpbpKxb+ISE2nn+UqOO3jHEdb7Cw+UyzL6qMCxxaHTIgXW+XUvdpvjEmxLKvEsiy1X/Ug5+8lqBkOV/sJsTGmwvkGYCo1k7N9gsb8i0hAc36EO1zrL3gHY0ypsztWlrMFq3iI8wrzssPuKB7TSIG/nJoLFuJ5693ctxFI9HSQtlDx73/cjfWvx5c+mhLpSM7hPinAKJujyKFyUYHjMc6fBf198A0l1CxOKJ7jqq3cDYnbhY/Np9SwHz/j/AiqAveTFyNRz2aRupYAo1Xw2MNZcLr6Zjf8fbXTuY+6znjGVCDMsqzTGtwfRs0wrHigQENNPMe5um+KhiPar05t5W5e5WHnWnobFf/+aQPuvxmj8JM2VSJt5Zw4l1K38LcsK1q9tD0qkppJvu4uVvR1/vewn2ZK2zXWp9y1+rWKfltU4L4DWVQj90vHKqTma++OT12g0LAf/5RPzWqADcUBeR7OIuJ1LMtKpKZjQ8Nf2FoHw4Ocb7SyG7myGQcU6lMZCWDufkdBTWcyrUbueXm4Xwgynpp5GD5Dxb8fcnZqiKwzO93VRm+XPj60VRg+Ni7QHzkn+EY5/x3tvMU5749p+tHSAdY734zVqtNHWyuSe4cwuwMEqEI3PxvJQKk6Mnmes37a5VrpF2prq2hfW+HXMsbYnUE6gHMs7VRqJgYBRPnaN6e/cA4vCQPGOP+7nJqhDDkaYuJ5lmU19Usv19lpRjzI+Qd0LDXj/PtSU/iP11V/+zh/b0VSc5W5lJohD1maf+FZdX42oObvR4kxJtu+ROL82ajli7WVin8RERERkQChYT8iIiIiIgFCxb+IiIiISIBQ8S8iIvL/7dzBbRNRFIXhcyUKsEQHLiG0EDoILYQOcAkodJB0EIUOSAvQAemAyB1cFh4kK0s2bzzv+zbWPG/O8l+8GYBJiH8AAJiE+AcAgEmIfwAAmIT4BwCASYh/AACYhPgHAIBJiH8AAJiE+AcAgEmIfwAAmIT4BwCASbwbPQAA3qqqL0mOST4mOSR5TXK7/P2+uw+DpgFcNPEPwKpU1V2Sr919rKqXJE9Jnrv7UFXXSX5U1WN3/xq7FODyuPYDwGr8i/vuPi5HL0mukvxenndJHoQ/wP+p7h69AQCSJFV1093fz55vk9x3dw2cBbAZ4h+A1aqqpyT77v4wegvAFrj2A8CaXSd5HD0CYCvEPwCrVFX7nO74P785vxoyCGADxD8Aq1FV90v0J8nnJDl/uXd5B+A4YBrAJoh/AFahqm5y+pb/bjn6s5zvz3533f0yZCDABnjhF4BVqKpdkrskP3OK/G/Lpz8/LWfp7odxCwEun/gHAIBJuPYDAACTEP8AADAJ8Q8AAJMQ/wAAMAnxDwAAkxD/AAAwCfEPAACTEP8AADAJ8Q8AAJMQ/wAAMIm/2tVF7wrSMaQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "The maximum error is: 0.049984407218554114\n", "The L2 error norm is: 0.004458860190283753\n" ] } ], "source": [ "# Step size and grid\n", "h = 0.1\n", "x = np.arange(0, 2*np.pi, h) \n", "\n", "# Target function\n", "y = np.cos(x) \n", "\n", "# Forward differences as an example - using diff \n", "# As always, please explore the documentation: https://numpy.org/doc/stable/reference/generated/numpy.diff.html\n", "y_forwarddiff = np.diff(y)/h \n", "\n", "# Consider where we compute the derivatives and translate grid\n", "x_forwarddiff = x[:-1:] \n", "\n", "# Reference(exact solution)\n", "y_exact = -np.sin(x_forwarddiff) \n", "\n", "# Plot results and look at the qualitative agreement\n", "plt.rcParams.update({'font.size': 22})\n", "plt.figure(figsize=(12, 8))\n", "plt.plot(x_forwarddiff, y_forwarddiff, \"ko\", markersize = 8, markerfacecolor=(0, 0.447, 0.741, 1), label = \"Forward difference approximation\")\n", "plt.plot(x_forwarddiff, y_exact, \"k-\", label = r\"Exact solution\")\n", "\n", "plt.legend(loc=\"upper left\")\n", "\n", "plt.xlabel(r\"$x$\")\n", "plt.ylabel(r\"$y$\")\n", "\n", "plt.grid()\n", "plt.show()\n", "\n", "# How well did we do (quantitavely)?\n", "errorLinfty = LA.norm(y_forwarddiff-y_exact, np.inf)\n", "print (\"The maximum error is: \", errorLinfty)\n", "\n", "errorL2 = LA.norm(y_forwarddiff-y_exact)/len(y_forwarddiff)\n", "print (\"The L2 error norm is: \", errorL2)\n", "\n", "# Consider why both of the above may be useful to learn about.\n", "# This may be a good opportunity to refresh your knowledge on norms and their mathematical foundations:\n", "# https://numpy.org/doc/stable/reference/generated/numpy.linalg.norm.html" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise 1**: using the code above, study the behaviour of the error for a variety of values of $h$. How would you best choose the latter? What do you infer from this investigation (and what do you expect)?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise 2**: Modify the above procedure developed for forward differences to accommodate for backward and/or central differences. Discuss its convergence behaviour. If you have time, attempt higher order derivative implementations as well" ] } ], "metadata": { "celltoolbar": "Slideshow", "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.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }