{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Forward Euler for single-variable model\n", "We will develop code to use the forward Euler method to solve a single variable ODE, then test it with the \"test\" equation (the degradation model).\n", "\n", "First, make a single-variable degradation model." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import scipy.integrate\n", "\n", "def degradation_model( t, y, params ):\n", " dydt = -params[0] * y;\n", " return dydt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Write the method\n", "\n", "The method should take in a function handle, the array of time steps, the initial value of y0 (a single state), and an array of params (though it will have only one entry).\n", "\n", "The outline is this\n", "\n", "- Make a 2D array to hold Y (1 x NT) \n", "- Set the value for the first time step\n", "- for each of the remaining time steps\n", " - compute dydt at current time step\n", " - set value of Y at next time step using Euler method\n", " \n", "The Euler method as we wrote in the notes is\n", "\\begin{equation*}\n", "\\widehat{y}(t + \\Delta t) = \\widehat{y}(t) + \\Delta t \\cdot f(t,\\widehat{y}(t))\n", "\\end{equation*}\n", "\n", "and since this is the code version, we really should add a vector of parameters to the function, i.e.\n", "\n", "\\begin{equation*}\n", "\\widehat{y}(t + \\Delta t) = \\widehat{y}(t) + \\Delta t \\cdot f(t,\\widehat{y}(t),\\vec{p})\n", "\\end{equation*}\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def forwardEuler1D( fn, t, y0, params ):\n", " \"\"\" [t,y] = forwardEuler1D( fn, t, y0, params )\n", " Given a function handle fn for dydt, find the value of\n", " y at timepoints t, given y(0) = y0.\n", " fn takes t, y, and params as input.\n", " t is 1-dimensional (NT,)\n", " y is 2-dimensional (1,NT) (NOTE WE ASSUME WE HAVE JUST ONE VARIABLE)\n", " \"\"\"\n", " # Write code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Compare the method's solution to the true solution\n", "\n", "Make an array of time steps starting at t=0, ending at t=4, with a time step of 0.1. Use y0 = 1 and a degradation rate constant of 2. \n", "\n", "Evaluate the true solution at those time steps, then use forward Euler to approximate the true solution at those time steps.\n", "\n", "Plot both solutions in an upper subplot and the difference between them in a lower subplot." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "# Write code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Evaluate the effect of step size on the correctness\n", "\n", "For values of $\\Delta t$ = 0.001, 0.01, 0.1, and 0.2, simulate the model and plot the difference between the forward Euler solution and the true solution." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "dts = np.array( [0.001,0.01,0.1,0.2])\n", "\n", "# Write code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### What happens if the step size is way too big?\n", "\n", "Use the forwardEuler method with $\\Delta t$ = 1. Plot the solution along with the true solution." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Write code here" ] }, { "cell_type": "code", "execution_count": 85, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEkCAYAAAA/7cqkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABjS0lEQVR4nO3dd3xO1x/A8c/JkAQRhBgJSWxZImLvVXtTq2ZRtVpdtFptaVFVq7Rq7xn92ZuEqJkSKwgiJBFEEiTIPr8/nngaZHsi67xfr+fV59577rnf56l8c3PvPd8jpJQoiqIouZ9edgegKIqi6IZK6IqiKHmESuiKoih5hEroiqIoeYRK6IqiKHmESuiKoih5RL5J6EKIZkKIXRncZ7AQomwG2lcTQpwUQkQLIb5IpZ2tEOK0EOKmEGKTEKJA4nqjxOWbidttkuzzdeL660KINhn5HIqi5A/5JqFn0mAg3QkdCAPGAbPSaPcLMEdKWQkIBz5MXP8hEJ64fk5iO4QQdkAfwB5oC/whhNDPQFyKouQDeTKhCyGmCCE+TbL8M1ADKCyEcBNCXBNCrBNCiMTtk4UQZ4UQl4UQi4VGT8AVWCeE8BZCmKR1XCnlQynlWSA2ldgE0AJwS1y1Cuia+L5L4jKJ21smtu8CbJRSRkspbwM3gTrp/DoURckn8mRCB5YDAwGEEHpozm4DgZrAp4AdUAFomNh+gZSytpTSATABOkop3QAvoL+U0llK+UIIMScxub/+mpiB2MyBx1LKuMTlQMAy8b0lEACQuP1JYnvt+mT2URRFAcAguwPIClJKfyFEqBCiJlAKOA+EAmeklIEAQghvwAY4DjQXQnwFFASKA1eAncn0O/6dfABFUZRMyJMJPdFSNNfAS6M5YweITrI9HjAQQhgDfwCuUsoAIcQPgHFyHQoh5gDNk9m0UUo5I51xhQJFhRAGiWfhVkBQ4rYgoBwQKIQwAMwS279c/1LSfRRFUYC8e8kF4H9obiDWBvan0u5l8n4khCgM9EyyLQIwfbkgpRyfePnl9Vd6kzlSUw3NPclxBgHbE9/vSFwmcfuRxPY7gD6JT8HYApWBM+k9pqIo+UOePUOXUsYIIdzRXK+OT7z/mVy7x0KIJcBl4D5wNsnmlcAiIcQLoL6U8kVqxxRClEZz3b0IkJB4Y9ZOSvlUCLEHGCalvAdMADYKIX5CczloWWIXy4A1QoibaJ6Y6ZMY4xUhxGbAB4gDRksp4zP2jSiKkteJvFo+N/Fm6Dmgl5TyRnbHoyiKktXy5CWXxOe2bwKHVTJXFCW/yLNn6IqiKPlNnjxDVxRFyY/ybEIXQowTQlxNHBHaVQgxOY32HkII1wz0byOE6Jdk2VEIsfIt4u0vhLgohLgkhDghhKiRQjtVB0ZRlGTl2YQOjAJaSyn7A1+hedZcl2wAbUKXUl4CrIQQ5TPZ322gqZTSEZgKLE6hnaoDoyhKsvJkQhdCLEIztH+vEGICEC2lfCSEMBVC3BZCGCa2K5J0GeglhDgjhPAVQjRObGMjhPAUQpxLfDVIbDsDaJw49P/lCNKdJD5qmFFSyhNSyvDExVNoBg+9/rlUHRhFUVKUJxO6lHIkcA/NqM6HaB5fREoZAXgAHRKb9gH+llK+LKZlIKWsg6bey/eJ6x6iOdN3AXoD8xPXTwQ8EwcWzUlc5wW8/EXQPIW6LyfS8RE+BPYms17VgVEUJUV5dmBREmWAkCTLS9FcgtkGDAGGJ9n2d+J//0VzSQXAEFgghHBGUy6gSirHekhiuV0ppTvgnNFghRDN0ST0RhndV1GU/C0/JPQXaGqiACCl/CfxMkozQF9KeTlJ25e1XuL577sZDzxAU35XD4hK5VjGicd7mZjnJNPmuZSygRBiNP/9MmkvpbwnhHBC8wunnZQyNJl9VR0YRVFSlCcvubzmKlDptXWrgfXAinTsbwYESykTgAHAyxuKr9R5SVQFTQkBpJTuKdR9aZC4fWGSdfcSb6b+DQyQUvomF4iqA6MoSmryQ0I/BtR8OZlFonVAMWBDOvb/AxgkhLgAVAOeJa6/CMQLIS4kuSnaHNidyTgno7nm/UfitXavlxuEEHvEf1PhTQA+S6z3Ys6rdWDME9d/huYaP1LKK8DLOjD7UHVgFCXPyhcjRYUQ84CdUspDics9gS5SygE6PIYRcBRolOSmpaIoyjuTH66hA0wD6gIIIX4H2gHtdXyM8sBElcwVRcku+eIMXVEUJT/ID9fQFUVR8gWV0NNJCFFRCFEyu+NQFEVJiUro6SCEaAic5M3HHxVFUXIMldDTIIRojWZU6QAp5clsDkdRFCVF+eUpl0wRQnQBlgDdpJTHszseRVGU1Kgz9BQk1jr/C80wfJXMFUXJ8VRCT4YQYgQwE2gppfw3u+NRFEVJD3XJ5TVCiM+BMUAzKeXN7I5HURQlvVRCT5RY6+V7oC/QREoZkMYuiqIoOYpK6GiT+SygFZpk/iCbQ1IURcmwfJ/QE+fX/BNwAppLKcOyOSRFUZRMydcJPXEu0ZVoZhlqnThFnaIoSq6UbxO6EMIY2Ihmirn2UsoX2RySoijKW8mXjy0KIQoBO4EYNIOGVDJXFCXXy3cJXQhRFDgABAB9pZQx2RuRoiiKbuSrhJ5YLfEI8C8wTE3FpihKXpJvEroQwhLNFHF7gU8SJ31WFEXJM/JFQk+c7f4YsFpKOUmqaZoURcmD8nxCF0JUR5PMf5NSzsjueBRFUbJKnn5sUQhRE9gDTJBSrs7ueBRFUbJStk0SXaJECWljY5Nl/UdGRnLr1i3Kly9PsWLFsuw4iqIo79K///77SEqZ7HSY2XaGbmNjg5eXV5b0ffjwYfr27cvOnTtp27ZtlhxDURQlOwgh7qS0Lc9dQ9+5cyd9+/bFzc1NJXNFUfKVXJ/QL1y4gLu7OwAbN25k+PDh7N69myZNmmRzZIqiKO9Wrr8pOnXqVNq0acOtW7f4/vvvOXjwII6OjtkdlqIoyjuXqxN6eHg4Bw8exMXFhcWLF+Ph4UHlypWzOyxFeUVsbCyBgYFERUVldyhKLmJsbIyVlRWGhobp3ifNhC6EWA50BB5KKR2S2S6AeUB74DkwWEp5Lt0RvIXNmzdjZWXF8uXL+frrr5k1axaff/45VapUeReHV5R0CQwMxNTUFBsbGzQ/LoqSOikloaGhBAYGYmtrm+790nMNfSWQ2t3FdkDlxNcINJNFvBM///wzN27c4P79+2zduhUHBwfKly//rg6vKOkSFRWFubm5SuZKugkhMDc3z/BfdWmeoUspjwkhbFJp0gXNkHoJnBJCFBVClJFSBmcokkxwdXWlU6dO9OjRgyJFimT14RQl01QyVzIqM/9mdPGUiyWaUrQvBSauyxKBNy9zcvFY4uPi+PvvvxkyZIhK5oqSitDQUJydnXF2dqZ06dJYWlpql2NidFM9eteuXdSsWZMaNWpgZ2fHX3/9lWr7H374gVmzZqXaxtvbmz179miXd+zYwYwZuqne8eLFC5o2bUp8fMoFVy9dusTgwYN1crx35Z3eFBVCjEBzWSbTl0YCT26m/r3VnJ99k6qjN1GwkErmipIac3NzvL29AU0iLVy4MF988YV2e1xcHAYGmU8FsbGxjBgxgjNnzmBlZUV0dDT+/v5vGbUmoXt5edG+fXsAOnfuTOfOnd+6X4Dly5fTvXt39PX1U2zj6OhIYGAgd+/ezTWXcnVxhh4ElEuybJW47g1SysVSSlcppWvJksmOXE1TvQFTOF31K2o8O0nAnFY8epDsoRRFScXgwYMZOXIkdevW5auvvnrjjNnBwUGblNeuXUudOnVwdnbmo48+euOsNiIigri4OMzNzQEwMjKiatWqAPj7+9OiRQucnJxo2bIld+/efSOWZs2aaUeNP3r0CBsbG2JiYpg8eTKbNm3C2dmZTZs2sXLlSsaMGZNqv4MHD2bcuHE0aNCAChUq4ObmluznX7duHV26dAFg4MCBbNu2Tbutf//+bN++HYBOnTqxcePGDH232UkXCX0HMFBo1AOeZPX187p9J3Gp4XysY/14sagld25czMrDKUqeFBgYyIkTJ5g9e3aKba5evcqmTZv4559/8Pb2Rl9fn3Xr1r3Spnjx4nTu3Blra2v69u3LunXrSEjQTDcwduxYBg0axMWLF+nfvz/jxo1LV2wFChRgypQp9O7dG29vb3r37v3K9tT6DQ4O5vjx4+zatYuJEye+0XdMTAx+fn68rCX14YcfsnLlSgCePHnCiRMn6NChA6C5T+fp6ZmumHOC9Dy2uAFoBpQQQgQC36OZWBkp5SI01QzbAzfRPLY4JKuCTarGewO5UbwsJXcNJmFde660W4l93Vbv4tCKkmk/7ryCz72nOu3TrmwRvu9kn+H9evXqleolB9DURfr333+pXbs2oLn2bGFh8Ua7pUuXcunSJQ4dOsSsWbM4ePAgK1eu5OTJk/z9998ADBgwgK+++irDcSYntX67du2Knp4ednZ2PHjw4I19Hz16RNGiRbXLTZs2ZdSoUYSEhLB161Z69OihvQRlYWHBvXv3dBLzu5Cep1z6prFdAqN1FlEGVHZtRXCxvci1Pamwpy+nw36jbruB2RGKouQ6hQoV0r43MDDQnlUD2sflpJQMGjSI6dOnp9mfo6Mjjo6ODBgwAFtbW+1Zb1qSHlsXg6+MjIy075OrJmtiYvLGcQYOHMjatWvZuHEjK1as0K6PiorCxMTkrWN6V3L1SFGAMhUdeTrKnaDFXXA9NY6joYE06f+1ekxMyZEycyb9LtjY2LBr1y4Azp07x+3btwFo2bIlXbp0Yfz48VhYWBAWFkZERATW1tbafSMjI/Hy8qJZs2aA5mbmy+0NGjRg48aNDBgwgHXr1tG4ceNkj/3vv/9Sp06dV655m5qaEhERkWy86ek3JcWKFSM+Pp6oqCiMjY0BzbX3OnXqULp0aezs7LRtfX19cXB4YzxljpXri3MBFClZlnLjD+NjWp+mN3/h6MJRxMbFZXdYipJr9OjRg7CwMOzt7VmwYIF2tLWdnR0//fQT7733Hk5OTrRu3Zrg4FdvkUkpmTlzJlWrVsXZ2Znvv/9ee3b++++/s2LFCpycnFizZg3z5s1749hffPEFf/75JzVr1uTRo0fa9c2bN8fHx0d7UzSp9PSbmvfee4/jx49rl0uVKkX16tUZMuTVK8bu7u7a6+m5QbZNcOHq6ip1XQ9dxsdyYfFHOD/YykmTZtiNWoeZaWGdHkNRMurq1atUr149u8NQkjh37hxz5sxhzZo1ADx//hxHR0fOnTuHmZkZANHR0TRt2pTjx4+/1WOdbyO5fztCiH+llK7Jtc8TZ+gvCX1DnEcu41K1T6n/wgP/uW0IzEU3NBRFeTdcXFxo3rw58fHxHDp0iOrVqzN27FhtMge4e/cuM2bMyLZknhl56gw9qRuHlmNz/AvuUobnvTbiaK9K6irZQ52hK5mVr8/Qk6rcaigPu2zAgjAsNnfi6NHD2R2SoihKlsqzCR3AsmYbEgbvQ19fn1pH+rFj6+pkH2NSFEXJC/J0Qgcws6mB6RgPnhhb0v7iJ2xeMoOYuIS0d1QURcll8nxCBzAqXo6y490JKuZK73sz2DZ3LGGR0dkdlqIoik7li4QOIIzNsB6zm7vluvB+5FpOzunLjXth2R2Worwz27ZtQwjBtWvXsi2Ge/fu0bNnT530tW3bNnx8fDK0z4MHD+jYsaO2zO/LSo4p8ff3T9fAomnTpr2y3KBBgwzFpSv5JqEDYFCA8kNXEew8jg7xh3nwV1c8L/tld1SK8k5s2LCBRo0asWHDBp31GZfBAXxly5ZNsQJiRmUmoU+ePJnWrVtz4cIFfHx8dFZf/fWEfuLECZ30m1H5K6EDCEGZrlN53Oo36otLFNvcjY1HzqibpUqeFhkZyfHjx1m2bNkr5WA9PDxo0qQJHTp0oGrVqowcOVJbV6Vw4cKMHz8ee3t7WrZsSUhICKApd/vpp5/i6urKvHnzOHz4MDVr1sTR0ZGhQ4cSHR3N2bNncXJyIioqimfPnmFvb8/ly5dfOeNduXIlXbt2pXXr1tjY2LBgwQJmz55NzZo1qVevHmFhmr+glyxZQu3atalRowY9evTg+fPnnDhxgh07dvDll1/i7OzMrVu3uHXrFm3btqVWrVo0btw42b9EgoODsbKy0i47OTkBmtGuX375JQ4ODjg6Or4xMvVlvC/L9wJ07NgRDw8PJk6cyIsXL3B2dqZ///7a7y61fj08PGjWrBk9e/akWrVq9O/fXzc5SEqZLa9atWrJ7Pbiyl754gcLGTC5gpyzfruMiYvP7pCUPMjHxye7Q5Br166VQ4cOlVJKWb9+fenl5SWllNLd3V0aGRnJW7duybi4ONmqVSu5ZcsWKaWUgFy7dq2UUsoff/xRjh49WkopZdOmTeXHH38spZTyxYsX0srKSl6/fl1KKeWAAQPknDlzpJRSTpo0SX7++edy1KhRctq0aVJKKW/fvi3t7e2llFKuWLFCVqxYUT59+lQ+fPhQFilSRP75559SSik//fRTbT+PHj3Sfo5JkybJ+fPnSymlHDRokDZWKaVs0aKF9PX1lVJKeerUKdm8efM3vod9+/ZJMzMz2axZM/nTTz/JoKAgKaWUbm5uslWrVjIuLk7ev39flitXTt67d++NeF9+B1JK2aFDB+nu7i6llLJQoUKvHOflckr9uru7yyJFisiAgAAZHx8v69WrJz09Pd+IN7l/O4CXTCGv5p4hUFnA2K4tCR/uo+iq7gy59hHT//yRcUMHU7RggewOTcmr9k6E+5d022dpR2iX+qWDDRs28MknnwDQp08fNmzYQK1atQCoU6cOFSpUAKBv374cP36cnj17oqenp61D/sEHH9C9e3dtfy/XX79+HVtbW23tl0GDBrFw4UI+/fRTJk+eTO3atTE2Nmb+/PnJxtW8eXNMTU0xNTXFzMyMTp06AZrKjRcvauY5uHz5Mt9++y2PHz8mMjKSNm3avNFPZGQkJ06coFevXtp10dFvPvjQpk0b/Pz82LdvH3v37qVmzZpcvnyZ48eP07dvX/T19SlVqhRNmzbV/pXxNlLqt0iRItSpU0f714KzszP+/v40atTorY6XrxM6gJ5VTQqP8uDpsq5MDPmaGfPu0+/Dz6hkoWrAKHlDWFgYR44c4dKlSwghiI+PRwjBr7/+Crw5GXFKlUqTrk9aejcloaGhREZGEhsbS1RUVLL7JC11q6enp13W09PTXp8fPHgw27Zto0aNGqxcuRIPD483+klISKBo0aLaqfZSU7x4cfr160e/fv3o2LEjx44dS3MfSLnEcGYl/ez6+voZvh+RnHyf0AEoZk2RUYeJWPU+kx/8xm8LH+DSZzLNq5fK7siUvCaNM+ms4ObmxoABA16ZuLlp06bamXjOnDnD7du3sba2ZtOmTYwYMQLQJEk3Nzf69OnD+vXrkz17rFq1Kv7+/ty8eZNKlSqxZs0amjZtCsBHH33E1KlTuX37NhMmTGDBggWZij8iIoIyZcoQGxvLunXrsLTUzEGftLxukSJFsLW1ZcuWLfTq1QspJRcvXqRGjRqv9HXkyBHq1atHwYIFiYiI4NatW5QvX57GjRvz119/MWjQIMLCwjh27Bi//vrrK0nbxsaGP/74g4SEBIKCgjhz5ox2m6GhIbGxsRgaGr5yvJT6zaonjfLfTdGUFCyO6bBdPK/cmc/FWgLWj2axh6+6Warkehs2bKBbt26vrOvRo4f2aZfatWszZswYqlevjq2trbZtoUKFOHPmDA4ODhw5coTJkye/0bexsTErVqygV69eODo6oqenx8iRI1m9ejWGhob069ePiRMncvbsWY4cOZKp+KdOnUrdunVp2LAh1apV067v06cPv/76KzVr1uTWrVusW7eOZcuWUaNGDezt7bXzgib177//4urqipOTE/Xr12fYsGHUrl2bbt264eTkRI0aNWjRogUzZ86kdOnSr+zbsGFDbG1tsbOzY9y4cbi4uGi3jRgxAicnJ+1N0ZfS068u5dniXJmWkEDs/m8xPL2QA/G1OGQ3jSk9a2NsmPpUXYqSkpxcnMvDw4NZs2ZpJ7dIqnDhwkRGRmZDVMpLqjjX29LTw7DdNGTbGbTWP0efq6MZvmg/D56+/dRYiqIoWUkl9BSIeh8j3l9NDYMAfno0no/nb8E74HF2h6UoOtWsWbNkz84BdXaeC6mEnhq7zugP3omlSTTL4r7hp79W8/e5wOyOSlEUJVkqoaelfF0Mhh2iiFlx1hn+xF63ZUzd5UNcvKrYqChKzqISenqUqIT+8EMYlrHnrwJziTn5FwOXnyHsWUx2R6YoiqKlEnp6FS6J3uBd6FVpw1TDlTQLWECX34/hc+9pdkemKIoCqISeMQUKQe+14PohI/R28l3MHHr/6cHOC2oiaiXnU+VzVflcAIQQbYUQ14UQN4UQE5PZbi2EOCyEuCiE8BBCWCXXT56gbwAdfoOW3/NewnE2Gs9k0gZPZuy9RnyCGoSk5FyqfK4qn4sQQh9YCLQD7IC+Qgi715rNAlZLKZ2AKcB0XQeaowgBjT+D7kuxi7/GQbNp7Dx6iiErz/L4ubquruQ8qnyuRr4vnwvUB/YnWf4a+Pq1NleAconvBfA0rX5zQvlcnfA7KuW0cvLZtIqy0zcLZaNfDsvLQY+zOyolB1Hlc1X53JxUPtcSCEiyHAjUfa3NBaA7MA/oBpgKIcyllKGZ/D2Te9g2gaH7KLiuJ/+L/4nxMZ/T489opnd3pFvNvHvlScmcX878wrUw3V7Drla8GhPqTEi1jSqfq6HK56bPF8ACIcRg4BgQBMS/3kgIMQIYAVC+fHkdHToHKGUHww6hv64X80Km8WfxcYzflMCFgCdM6lAdQ31171nJPqp87qvye/ncIKBckmWrxHVaUsp7aM7QEUIUBnpIKR+/3pGUcjGwGDTFuTIXcg5VpCwM2YvYPIBRfrOxrzCMQSeac+XeExb2c8GiiHF2R6jkAGmdSWcFVT73P6p8LpwFKgshbIUQBYA+wI6kDYQQJYQQL/v6Gliu2zBzCeMi0G8L1OhL03tLOVp1K9eCwuj4+3H+vROW3dEp+ZQqn/sfVT4XEEK0B+YC+sByKeXPQogpaC7O7xBC9ETzZItEc8lltJTyzQtYSeTY8rm6ICW4T4NjM4ks15SeoSO5+Ri+62jHwPrWKf5Jq+RNqnyuklkZLZ+brmvoUso9wJ7X1k1O8t4N0M3DpXmBENBiEphZUXjXeHZbhPFF8W/5fscVzt0NZ3p3RwoWUJNFKYqiW+puXVaqNQj6bUI/zI/ZEZ/zc0MDdly4R5cF/3DzoTrzUbKfKp+bt6iEntUqt4YhuxFxMfS/MpztHSD0WQxdFhxn98Xg7I5OUZQ8RCX0d6FsTRh2CAqXwsl9MIfee0SV0qaMXn+Oqbt8iFWleBVF0QGV0N+VYtYwdD9YulJ870i2OHoxuL41y47fpu/iU2qKO0VR3ppK6O9SweIw4H9g1xWDw5P5wXAV8/s44RP8lA7zPTlx81F2R6goSi6mEvq7ZmgMPVdA/TFwZjGdr01kxwgXzEwM+WDZaeYfvkGCqtqoZAFVPleVz1Wygp4etPkZ2v4C13ZTaV8/dgytRucaZZl90JdBK84QGpnqY/yKkmGqfK4qn6tkpXoj4f3VcP8Shda0Y07rIkzv7sjp22F0mH+cs/5qdKmiG6p8rka+L5+bVa88Uz5XF+6cknKGjZS/2EoZcFZeDnosm8w8Iit8vVsu8rgpExISsjtC5S2o8rmqfG5OKp+rZLXydeHDg7CuB6zsiH3PZewc+x4T3C4yfe81zvqHMatXDYoWLJDdkSpv6f60aURf1e01bKPq1Sj9zTeptlHlczVU+Vzl3ShRCT48BOvfh00fUKTdTP7oP4xVJ/z5ec9VOsw/zu/9auJSvlh2R6rkMqp87qvye/lc5V0pXBIG7wK3D2HPF4gngQxu+T3O5YsxZv053l90kq/aVmVYowro6akCX7lRWmfSWUGVz/2PKp+rvFsFCkHvteD6IfwzF/4ejnMZE3aPa0xru1JM23OND1edJeyZmrtUSR9VPvc/qnxuFsnT5XN1QUo4PgcO/wg2jaH3WqSxGWtP3WHqrqsUK2TI/D41qVvBPLsjVdKgyucqmZXR8rnqDD2nEgIafwbdl8LdU7C8DeJJIAPq2/C/0Q0oWMCAvktO8fvhG8SrgUiKoqASes7n1AsG/A1P78HSVhB8EfuyZuwc24jONcry20FfBi4/zcMIVQtGyThVPjdvUQk9N7BtoinspacPK9rBzcMUNjJgTm9nZvZ04t874bSf54nH9YfZHamiKNlIJfTcopSdpgRvMRvNo43n1yGE4H3Xcuwc04gShY0YvOIsP+/2ISZOleNVlPxIJfTcpEhZGLIXbBrB9lHgMQOkpHIpU7aNbsjA+tYs8bxNjz9PcPvRs+yOVlGUd0wl9NzGuAj02wI1+oLHdNgxBuJjMTbUZ0oXB/4aUIu7Yc/pMN+Trf8GZne0iqK8Qyqh50YGBaDrn9DkKzi/Ftb3hmjNAIs29qXZ+0ljHCzN+HzLBcZv8iYiKjabA1ZyAlU+V5XPVXIqIaDFJOg0H/w8YEV7iLgPQNmiJmwYXo/PWldhu3cQHX8/zvm74dkbr5LtVPlcVT5XyelqDYJ+myD0FixtDSHXAdDXE4xrWZlNH9UnLl7Sc9FJFhxRz6znV6p8roYqn6vK5+YOQeeknFlJyunlpLx9/JVNj5/HyNHr/pXWE3bJXn+ekIHhz7MpyPxJlc9V5XNzVPlcIURbYB6gDyyVUs54bXt5YBVQNLHNRCnlnrf/daOkW9mamsca1/WENV2h2yJw6AGAmYkhv/etSfOqFkzefpm2c48xrZsjnWqUzd6Y8yHPzb48CtDtgJ0S5QrT+P0qqbZR5XM18n35XCGEPrAQaA0EAmeFEDuklEkvXn0LbJZS/imEsAP2ADZvFZmSccWsNQOQNvYHt6HwJAgajAUhEELQo5YVrjbF+HSTN2M3nMf9+kOmdHGgsJEqupmXqfK5r8rv5XPrADellH4AQoiNQBcgaUKXQJHE92bAvbeOTMmcgsVhwP/gfx/Bwe/gSQC0naEZZQpYmxdi80f1+f3wDRa438TLP5y5fZxVnfV3JK0z6aygyuf+R5XPBUsgIMlyYOK6pH4APhBCBKI5Ox+rk+iUzDE0hp4roP4Y7v+7jKnrWhAR+V9ZAEN9PT57ryobR9QnPkHSa9FJ5hz0JTZejTDNi1T53P/k+/K5QoieQFsp5bDE5QFAXSnlmCRtPkvs6zchRH1gGeAgpUx4ra8RwAiA8uXL17pz545OP4zypp37x/Nt8EFKF7RgetNZuJRyeWX706hYvt9+hf+dD6JGuaLM7e2MbYm0/5xW0k+Vz1UyKyvK5wYB5ZIsWyWuS+pDYDOAlPIkYAyUeL0jKeViKaWrlNK1ZMmS6Ti0kl4XLlxgx44dLFq0iAkTJlC3bl2uXbtGpzZzWNXyD/T0jRiyfwjzz80nNv6/gUZFjA2Z09uZBf1q4v/oGe3nebLu9B3dPEKlKMo7lZ4zdAPAF2iJJpGfBfpJKa8kabMX2CSlXCmEqA4cBixlKp2rCS5059GjR3Tp0oXixYsTGxvL0aNHEUJw584dXv7ifBb7jJlnZ/L3jb+xM7djeuPpVDCr8Eo/959E8aXbBTxvPKJlNQtm9HCipKlRcodUMiAnn6ErOZvOz9CllHHAGGA/cBXN0yxXhBBThBCdE5t9DgwXQlwANgCDU0vmim6VKFECNzc3jI2NuXnzJlZWVixbtoykfwUVMizEjw1+ZE6zOQRFBtF7Z282X9/8ypl4aTNjVg2pw/ed7PC8+Yi2c49x0OdBdnwkRVEyIV0jRaWUe6SUVaSUFaWUPyeumyyl3JH43kdK2VBKWUNK6SylPJCVQSv/SUhI4K+//sLJyYmqVasyatQoypUrR58+fZJt38q6FX93/huXUi5MPTWVMUfG8OjFI+12PT3BkIa27BrbiFJFjBm+2osJbhdVPZi3pM5vlIzKzL8ZNadoLubj48OIESOIj49nyZIlODg4YGFhgYeHB3Z2dqnumyAT2HBtA7O9ZlO4QGF+bPAjzco1e6VNTFwCcw/5sujoLcqYmTCrVw3qV1RzmGbU7du3MTU1xdzcPMVnvBUlKSkloaGhREREYGtr+8q21C65qISeC0VFRTF9+nT++OMPfvjhBz7++GP09DR/bIWFhVG8ePF093Uz/CYTPSdyPfw6var04gvXLyhoWPCVNv/eCefzzd74hz5naENbvmpbFWNDfZ1+prwsNjaWwMDAtx6IouQvxsbGWFlZvfFsu0roecjRo0f56KOPsLe3Z/78+dpBFm8jJj6GBecXsPLKSqyLWDOj8QzsS9i/0uZ5TBwz9l5j9ck7VCxZiDm9nXGyKvrWx1YUJWPe9rFFJQcICwtj2LBhfPDBB8yYMYOtW7fqJJkDFNAvwGeun7H0vaW8iHvBB3s+YPHFxcQnxGvbFCxgwJQuDqz5sA7PouPp9scJNRhJUXIYldBzOCklGzduxN7eHhMTE65cuULXrl2z5Fh1ytRha+ettLJuxe/nf2fI/iEERrw661HjyiXZ/2kTOtcoy7zDN+j+xwmu34/IkngURckYdcklB3v+/Dk9e/YkICCAJUuWUK9evXdyXCklu2/v5udTPyORfFP3GzpV6PTGDb29l4KZtO0ykVFxfNq6MiMaV8BAX50jKEpWUpdccqmEhAR69uzJuXPn3lkyB01VvY4VOrK181aqFqvKpOOT+OLoFzyJfvJKu3aOZTgwvgktq1swc991eiw6yc2H6mxdUbKLOkNXUhWfEM+KKytYeH4hxU2K83Ojn6lX5tVfLlJKdl0M5rvtl3keE8+X71VlaCNb9PXUI3qKomvqDF3JNH09fYY5DmNdh3UUMizE8APD+fXsr0TH/zd5gBCCTjXKcmB8E5pVKcnPe67y/l8n8QtRhZ0U5V3KUwndw8ODjh07ZmiflStXcu9e+su3SykZN24clSpVwsnJiXPnziXbbsOGDTg6OuLk5ETbtm159EgzGjMsLIzWrVtTuXJlWrduTXh4eIb6zS525nZs6riJ3lV7s9pnNX1398U33PeVNhamxvw1oBbz+jhz82Ek7eZ5stTTT81jqijvSJ5K6JmR0YS+d+9ebty4wY0bN1i8eDEff/zxG23i4uL45JNPcHd35+LFizg5OWmL+8+YMYOWLVty48YNWrZsqZ11PD39ZjcTAxO+rfctC1suJPRFKH139WWNzxoSklRJFkLQxdmSg+Ob0LhyCX7afZVei05w86E6W1eUrJZrE/rkyZOZO3eudnnSpElcuHCByMjIZGfSnjJlCrVr18bBwYERI0YgpcTNzQ0vLy/69++Ps7MzL168SPO427dvZ+DAgQghqFevHo8fPyY4OPiVNi8nbH327BlSSp4+fUrZsmW1+w8aNAjQzL+4bdu2dPebUzSxasLfnf+mQdkGzDw7k48OfsSDZ68W8bIoYsySga7M7e2M36NntJ/vyZ8et4hTz60rStZJafborH7VqlXrjdmsM+L27duyZs2aUkop4+PjZYUKFaSbm1uKM2mHhoZq9/3ggw/kjh07pJSaGczPnj2r3fbpp5/KGjVqvPGaPn26lFIz03fS2blbtGjxyv4vbdmyRZqamsrSpUvLxo0by7i4OCmllGZmZto2CQkJ2uX09puTJCQkyM3XN8vaa2vLBusbyP239yfb7uHTKDlyjZe0nrBLdvrdU14NfvKOI1WUvAPwkink1Vx7hm5jY4O5uTnnz5/nwIED1KxZE3Nzc+1M2np6etqZtAHc3d2pW7cujo6OHDlyhCtXriTb75w5c/D29n7jNXHixHTHFhsby59//sn58+e5d+8eTk5OTJ8+/Y12InHy5txKCEGvKr3Y3HEz5U3L8/nRz/n2+LdExrx6eaWkqRF/flCLhf1cCAp/QaffjzPv0A1i4tTZuqLoUq5N6ADDhg1j5cqVrFixgqFDhwLJz6QdFRXFqFGjcHNz49KlSwwfPjzFQknjx4/H2dn5jdfLa92WlpYEBPw3xWpgYOAbQ/BfzjxesWJFhBC8//77nDhxAoBSpUppL6UEBwdjYWGR7n5zKhszG1a3X81HTh+x028nPXf25PzD82+06+BUhoOfNaWdQxnmHPKl84LjXAp8kkyPiqJkRq5O6N26dWPfvn2cPXuWNm3apNjuZfIuUaIEkZGRuLm5abclnTkc0j5D79y5M6tXr0ZKyalTpzAzM6NMmTKvHM/S0hIfHx9CQkIAOHjwoHbWkc6dO7Nq1SoAVq1aRZcuXdLdb05mqGfImJpjWNVW89kG7xvM7+d/Jzbh1TrqxQsVYH7fmiweUIuwZzF0WXic6Xuu8iImPrluFUXJAIPsDuBtFChQgObNm1O0aFH09VMu51q0aFGGDx+Og4MDpUuXpnbt2tptgwcPZuTIkZiYmHDy5ElMTExSPWb79u3Zs2cPlSpVomDBgqxYsUK7zdnZGW9vb8qWLcv3339PkyZNMDQ0xNrampUrVwIwceJE3n//fZYtW4a1tTWbN29Os9/cxNnCGbdObsw4M4PFFxdzIugE0xtPx8bM5pV279mXpm4Fc6bvucpfx/zYd+U+07s70qDiG1PRKoqSTrl6pGhCQgIuLi5s2bKFypUr6ygyRVcO+B/gx5M/EpsQy5e1v6Rn5Z7J3jM4cesRX/99iTuhz+lbpxwT21XHzMQwmR4VRcmTI0V9fHyoVKkSLVu2VMk8h3rP5j3+7vw3NUrWYMrJKYw7Mo7QF6FvtGtQsQT7PmnCiCYV2HQ2gPfmHOXAlfvZELGi5G65+gxdyR0SZALrrq5j7r9zKVygMFMbTqWJVZNk214MfMxXbhe5dj+CDo5l+L6THRZFjN9xxIqSc+XJM3Ql99ATegywG8CGjhswNzFn9OHR/HTqJ17EvTmQy8mqKDvHNuLLNlU5ePUBLWcfZf3puySo8gGKkqZcndDnz59P9erV6d+/P9u2bWPKlCmptm/WrBkZ+avA39+f9evXa5cvXbrE4MGDMxsu165do379+hgZGTFr1qwU20kpmTRpElWqVKF69erMnz9fuz6lei+rVq2icuXKVK5cWfsUTU5TpVgVNnTYwCC7QWy6von3d77PldA3xwMY6usxunkl9n3SGIeyZnzzv0v0XqxK8ypKmlIacZTVr7cdKSqllFWrVpUBAQFSSinr168vQ0JCUm3/+qjQtLi7u8sOHTq8sq5ly5byzp07GQ9WSvngwQN55swZ+c0338hff/01xXbLly+XAwYMkPHx8dr9pJRy9+7dsm3btjIhIUGePHlS1qlTR0qpGQVra2srQ0NDZVhYmLS1tZVhYWGZivFdOXnvpGyxuYV0XuUsl1xcIuPi45Jtl5CQIDefvStr/LhfVvpmt5x94LqMik2+raLkB+TFkaIjR47Ez8+Pdu3a8csvv2BkZESJEiWIiIjA1taW2FjN889Pnz59ZXnLli3UqVOHKlWq4OnpCWjOxBs3boyLiwsuLi7aQUATJ07E09MTZ2dn5syZA0CnTp3YuHFjpmK2sLCgdu3ab8zi/bo///yTyZMno6enp90PUq73sn//flq3bk3x4sUpVqwYrVu3Zt++fZmK8V2pV6Yef3f+mxblWzDv3DyG7h9KUGTQG+2EEPRyLcehz5rSwbEM8w7foN08T077vXlzVVHyu1yb0BctWkTZsmVxd3fHwsICFxcXQDNQqFmzZuzevRuAjRs30r17d20SjYuL48yZM8ydO5cff/wR0CTMgwcPcu7cOTZt2sS4ceMATWXExo0b4+3tzfjx4wFwdXXV/iJwd3dPdlRpgwYN3uqz3bp1i02bNuHq6kq7du24ceMGAEFBQZQrV07bzsrKiqCgoBTX53RmRmbMajqLnxv9zPXw6/Tc0ZOdt3ZqC6olVaKwEXP71GTV0DrExCXQe/EpvnK7QPizmGyIXFFypnQNLBJCtAXmAfrAUinljNe2zwGaJy4WBCyklEV1GGeqgoODKVmypHZ52LBhzJw5k65du7JixQqWLFmi3da9e3cAatWqpa3zEhsby5gxY/D29kZfXx9f31frfCdlYWGhLbfbvHlz7TB/XYqOjsbY2BgvLy/+/vtvhg4dqv0lktcIIehcsTMuFi58c/wbvjn+DZ6BnkyqNwkzI7M32jetUpID45sw79ANlh6/zaGrD/mmfXV6uFjm6ro4iqILaZ6hCyH0gYVAO8AO6CuEsEvaRko5XkrpLKV0Bn4H/s6CWFNkYmLySm2Whg0b4u/vj4eHB/Hx8Tg4OGi3vaz18rLOC2iG+5cqVYoLFy7g5eVFTEzKZ31RUVHa0aRpnaEvXLhQuy4jNdetrKy0v3i6devGxYsXgZTrveTmOjAvWZlasaLNCsbVHMfBOwfpsaMHZ4LPJNu2YAEDvm5fnV1jG2FjXpAvtlygz+JTqua6ku+l55JLHeCmlNJPShkDbAS6pNK+L7BBF8GlV/Xq1bl58+Yr6wYOHEi/fv0YMmRImvs/efKEMmXKoKenx5o1a4iP19QVeb3OC4Cvr6/2F8TLM/TXXy+vwY8ePVq77mU99PTo2rUr7u7uABw9epQqVaoAKdd7adOmDQcOHCA8PJzw8HAOHDiQam2bnEpfT5/hTsNZ234tJgYmDDswjN+8fiMmPvlfsNXLFMFtZAOmdXPkavBT2s07xm8HrhMVq+rCKPlUSndLX76Anmgus7xcHgAsSKGtNRAM6KfVry6ecrG2tpYhISHy2bNn0s7OTiYkJGi3BQcHS2NjYxkeHq5dl/Qpl5CQEGltbS2llNLX11c6OjpKJycn+dVXX8lChQpJKaWMiYmRzZs3l05OTnL27NlSSilHjx6traWeUcHBwdLS0lKamppKMzMzaWlpKZ880dQGb9eunQwKCpJSShkeHi7bt28vHRwcZL169aS3t7eUUvPEx6hRo2SFChWkg4PDK0/sLFu2TFasWFFWrFhRLl++PFPx5STPYp7JKSemSIeVDrLH9h7yRtiNVNs/fBolP9lwTlpP2CWbzDwij/k+fEeRKsq7RSpPuaQ5UlQI0RNoK6Uclrg8AKgrpRyTTNsJgJWUcmwKfY0ARgCUL1++1p07dzLzOyhZn3zyCZ06daJVq1YAuLm5sX37dtasWaOzY0RHR9O0aVOOHz+OgUGurmuWaxwNOMrkE5OJjInkM9fP6FutL3oi5T8sj994xHfbL3P70TM6OpXh2w52lDZTI02VvCO1kaLpSej1gR+klG0Sl78GkFK+MWODEOI8MFpKeSKtoHQ99P/BgwecPn2azp07M3bsWPbu3cuePXu0lyt04caNGwQFBdGsWTOd9amk7dGLR3x/4nuOBR6jYdmGTGk4BYuCFim2j4qN56+jfiz0uImhnmB86yoMbmCDgX6ufahLUbTeNqEbAL5ASyAIOAv0k1Jeea1dNWAfYCvT6hRVy0XJGCklW3y38OvZXzE2MOb7+t/TyrpVqvvcCX3GDzuu4H49hGqlTfmpqwOuNsXfUcSKkjXeqpaLlDIOGAPsB64Cm6WUV4QQU4QQnZM07QNsTE8yV5SMEkLwftX32dxpM2ULl2W8x3gm/zOZZ7HPUtzH2rwQywfXZtEHtXj6Ipaei07y5ZYLhEZGv8PIFeXdUdUWlVwnNj6WPy/8ybLLyyhbqCzTG0/H2cI51X2ex8Qx//BNlnr6UcjIgC/bVKVvnfLo66ln15XcRVVbVPIUQ31DxrmMY0WbFUgkg/YNYqH3wjemu0uqYAEDJrarxt5PGlO9jCnfbrtMl4XH+fdO+DuMXFGylkroSq7lUsqFLZ220LFCRxZdWMTgvYO5+/RuqvtULmXKhuH1mN+3JiER0fT48wRfbrnAI3UZRskD1CUXJU/Y57+PqSenEpsQy4TaE+heuXuapQCeRccx/8gNlh+/jbGhPp+1rsKAetbqaRglR3urp1yyikroiq7df3afb49/y+n7p2lRrgU/NPiBYsbF0tzv5sNIftx5Bc8bj6hW2pQfO9tTt4L5O4hYUTJOXUNX8oXShUqz+L3FfOH6BZ5BnnTf0Z3jQcfT3K+SRWFWD63Dog9ciIiKo/fiU4zbcJ7gJ2/OqKQoOZk6Q1fypOth15noOZGbj2/St1pfPqv1GcYGaY8YfRETz58eN1l0zA99IRjdvCLDGlfA2FD/HUStKGlTl1yUfCk6Ppq5/85l7dW1VDCrwIzGM6huXj1d+waEPefn3VfZd+U+5Yqb8G0HO96zK6VK9CrZTl1yUfIlI30jJtSZwF+t/yIyJpJ+e/qx/PJy4hPSrsZYrnhBFg2oxbphdTEx1OejNf8yYNkZbjxQ85oqOZc6Q1fyhcdRj5lyagoH7xzEtZQr0xpNo0zhMunaNzY+gbWn7jDnoC/PYuIZWN+aT1tWwaxg6lMJKkpWUJdcFAVNPZjtt7Yz/fR09IU+k+pNokOFDunePzQymlkHfNl49i5mJoZ81roK/eqUV485Ku+USuiKkkRARADfeH6Dd4g37W3bM6neJIoUKJLu/a/ce8LUXT6c8gujSqnCfNfRjsaVS6a9o6LogEroivKauIQ4ll5ayqILiyhZsCTTGk2jduna6d5fSsn+Kw+Ytucqd8Oe06q6Bd+0r06FkoWzMGpFUQldUVJ0KeQSEz0nEhARwBCHIYxxHoOhfvqvjUfHxbPiH38WHLlJdFw8g+rbMLZlZcxM1PV1JWuohK4oqXge+5xfvX7FzdeNasWrMd3pGypZ18xQHw8joph9wJdNXgEUNTHk01ZV6Fe3PIbq+rqiY+qxRUVJRUHDgnxf/3vmN59PYZ9AIjv048K0iSQ8S7nW+ussTI2Z0cOJXWMbUb1MEb7fcYU2c49xyOcBaooA5V1RCV1RgNjYWJqXb86MXkt5WMuaAqu3c7NtWx5v3YqMT/u59Zfsy5qxblhdlg7UnEANW+1F/6WnuXLvSVaFriha6pKLku/ExsYyduxY/Pz8uHfvHnfu3CEyMpLVq1czYMAAAJ6fP8/DGb/w4sIFjKpVo9TECRSqVy9jx4lPYP3pu8w95MvjF7H0dLHiizZVKVVETVqtZJ66hq4oSSQkJPD3338TExPD8uXLuXDhAhYWFly6dAk9vf/+aJVSErF3Lw9n/UbsvXsUbt4ciy+/xKiCbYaO9+RFLAvdb7Lin9sY6OkxvLEtI5pWpLCRga4/mpIPqGvoivKakJAQPvnkE2rWrImenh5r1659JZmDZh7TIu3bU2HvHkp+9hnPz5zBr3Nn7v88jbjw9M90ZGZiyDftq3Pos6a0qG7B/CM3afarB2tP3SEuPkHXH03Jx9QZupKvXLlyhREjRgCwePFifv/9dwwNDfn999/T3Dfu0SNCfl/A4y1b0CtcmBKjPqZ4v36IAgUyFMP5u+FM23OVs/7hVCxZiIntqtOquoUq/KWkizpDV/K9qKgovvvuO5o1a8YHH3yAp6cn9vb2GBgYMHXq1HT1YVCiBGV+/AHbbf/DxNGRhzN+4VanTkQcOpShJ1lqli/G5o/q89eAWkgJw1d70XvxKbwDHmfy0ymKhjpDV/I8d3d3PvroIxwdHZk/fz6WlpZv3aeUkmeenjz4ZSYxt25RsHZtLCZOwMTePkP9xMYnsPFsAPMO+fIoMoYOTmX48r2q2JQo9NYxKnmTuimq5EuhoaF8+eWXHDx4kAULFtClSxedH0PGxfF4yxZC5v9O/OPHmHXpQsnxn2JYqlSG+omMjuOvo7dY6nmb2PgE+tUtz9gWlSlpaqTzmJXcTV1yUfIVKSXr16/HwcGBwoULc+XKlSxJ5gDCwIBifftS8cB+zD8cytPdu7nVth0hvy8g4fnzdPdT2MiAz9+rytEvm9G7djnWnb5L01/dmXPQl8jouCyJXcl70nWGLoRoC8wD9IGlUsoZybR5H/gBkMAFKWW/1PpUZ+hKVrh9+zYff/wx9+7dY8mSJdStW/edHj8mMJCHv/1GxN59GFhYUHL8eMy6dEboZezcyS8kklkHrrPn0n3MCxVgXMvK9K1TngIG6hwsv3urSy5CCH3AF2gNBAJngb5SSp8kbSoDm4EWUspwIYSFlPJhav2qhK7oWnR0NDY2NnzyySd8/vnnGBpmX4Gs5+fO8WD6DKIuXcLYzg6LiRMoVKdOhvs5fzecGXuvcfp2GOWLF+Tz96rQyaksenrqiZj86m0Ten3gByllm8TlrwGklNOTtJkJ+Eopl6Y3KJXQlawQHx+Pvn7OmNBZJiTwdPceHs6eTVxwMKatW2Hx+ecUsLHJWD9S4uEbwi97r3HtfgTVSpvyVduqNK+qHnXMj972GrolEJBkOTBxXVJVgCpCiH+EEKcSL9EoyjuXU5I5gNDTw6xTRyru3UPJTz8h8p8T3OrUmQfTZxD/JP21XYQQNK9qwZ5xjZnXx5kXsfEMXelFr0UnOXM7LAs/gZLb6OqCnAFQGWgG9AWWCCGKvt5ICDFCCOElhPAKCQnR0aEVJWfTMzamxMiRVNq/j6JduxC2ejW33mtD2Oo1yNjY9PejJ+jibMmhz5ryU1cH7oY95/2/TjJkxRlV/EsB0pfQg4BySZatEtclFQjskFLGSilvo7nmXvn1jqSUi6WUrlJK15Il1ZRd2c3Dw4OOHTtmaJ+VK1dy7969dLeXUjJu3DgqVaqEk5MT586dS7ZdTEwMI0aMoEqVKlSrVo2tW7e+sn3r1q0IIUh6mW769OlUqlSJqlWrsn///gx9juxgULIkZaZOxfZ/f2NkV50H06bh16kzEUfcMzQwyVBfjw/qWXP0y+ZMbFeNc3cf02H+ccZuOI9fSGQWfgIlp0tPQj8LVBZC2AohCgB9gB2vtdmG5uwcIUQJNJdg/HQXppJTZDSh7927lxs3bnDjxg0WL17Mxx9/nGy7n3/+GQsLC3x9ffHx8aFp06babREREcybN++VJ1Z8fHzYuHEjV65cYd++fYwaNYr4DJS5zU7G1apRfvlyrBb9CUIQOGoUd4cMJerq1Qz1Y1JAn5FNK3Lsq+aMbl6RQz4PaD3nGF+5XSAwPP2PTCp5R5oJXUoZB4wB9gNXgc1SyitCiClCiM6JzfYDoUIIH8Ad+FJKGZpVQSsZN3nyZObOnatdnjRpEhcuXCAyMpKePXtSrVo1+vfvrz1TnDJlCrVr18bBwYERI0YgpcTNzQ0vLy/69++Ps7MzL168SPO427dvZ+DAgQghqFevHo8fPyY4OPiNdsuXL+frr78GQE9PjxIlSmi3fffdd0yYMAFjY+NX+u3Tpw9GRkbY2tpSqVIlzpw5k9mv550TQmDarBkVdmyn1LffEn3tGre79+DepEnEPkz1AbE3mJkY8mWbanhOaM6g+jZs875H81keTN5+mQdPo7LoEyg5UbquoUsp90gpq0gpK0opf05cN1lKuSPxvZRSfialtJNSOkopN2Zl0ErGDR06lNWrVwOa8rEbN27EysqK8+fPM3fuXHx8fPDz8+Off/4BYMyYMZw9e5bLly/z4sULdu3aRc+ePXF1dWXdunV4e3tjYmLC+PHjcXZ2fuM1Y4ZmqEJQUBDlyv13xc7KyoqgoFev2D1+/BjQJG4XFxd69erFgwcPADh37hwBAQF06NDhlX3S029uIAwNKf5Bfyru30fxwYN5smMnt9q249Gff5KQjl+YSZUobMTkTnYc/bIZvVzLsf70XZrMdGfanquEPYvJok+g5CRqlEI+YWNjg7m5OefPn+fAgQPUrFkTc3Nz6tSpg5WVFXp6ejg7O+Pv7w9o6p/UrVsXR0dHjhw5wpUrV5Ltd86cOXh7e7/xmjhxYrpji4uLIzAwkAYNGnDu3Dnq16/PF198QUJCAp999hm//fabLr6CHE3fzIxSE76i4u5dFG7UiJB587nVrj1PduxAJmSsxG4ZMxOmdXPk8OdN6eBYhiWefjT+5Qi/HbjOk+fpvwmr5D4qoecjw4YNY+XKlaxYsYKhQ4cCYGT0X60QfX194uLiiIqKYtSoUbi5uXHp0iWGDx9OVFTyf7qndYZuaWlJQMB/T70GBga+URzL3NycggUL0r17dwB69erFuXPniIiI4PLlyzRr1gwbGxtOnTpF586d8fLySle/uVGB8uWxmj8P6zWrMTA3595XE/B/vzfPMzFmw9q8ELN7O3Pg0yY0rVqS34/cpNEvR5h7yJenUSqx50lSymx51apVSyrvVnR0tKxSpYq0tbWVcXFx0t3dXXbo0EG7ffTo0XLFihUyPDxcWlhYyOfPn8uIiAhpb28vv//+eymllB07dpRHjhxJ9zF37dol27ZtKxMSEuTJkydl7dq1k23Xu3dvefjwYSmllCtWrJA9e/Z8o03Tpk3l2bNnpZRSXr58WTo5OcmoqCjp5+en/Ux5SUJ8vHy8bZv0bdJU+lStJgPGjpPRd+9muj+fe0/kiNVnpfWEXdLx+33y98O+MiIqVocRK+8C4CVTyKtqDqx8pECBAjRv3pyiRYumOgCnaNGiDB8+HAcHB0qXLk3t2rW12wYPHszIkSMxMTHh5MmTmJiYpHrM9u3bs2fPHipVqkTBggVZsWKFdpuzszPe3t4A/PLLLwwYMIBPP/2UkiVLvtIuOfb29rz//vvY2dlhYGDAwoULc9SgIl0QenqYdemC6XvvEbpiBaFLlhLp7k6xAQMoMfIj9IsUyVB/1csU4a8BrlwOesLcQ77MOuDLsuO3GdGkIgPrW1NITYmX66nyuflIQkICLi4ubNmyhcqV3xgmoORwsQ8eEjJvHk/+9z/0zcwoMXYMxXr3RhhkLhFfCHjMnEO+eFwPwbxQAYY3qcCAeiqx53SqfK6Cj48PlSpVomXLliqZ51KGpSwoO+1nbLe6YVSlCg+m/oRfl65EHj2aoYFJL9UoV5SVQ+rw96gG2JUtwoy912g8050/PW7xTJXszZXUGbqi5EJSSiLd3Xn4y0xi7tyhUIMGWEyYgHHVKpnu89zdcOYdusFR3xCKFTRkeJMKDKxvQ2F1xp6jqDN0BYD58+dTvXp1+vfvz7Zt25gyZUqq7Zs1a0ZGfun6+/uzfv167fKlS5cYPHhwZsPl2rVr1K9fHyMjI2bNmpViu/79+1O1alUcHBwYOnQosYn1UZ48eUKnTp2oUaMG9vb2r1yXX7VqFZUrV6Zy5cqsWrUq0zFmFyEEpi1aUGHnDkp98w0vrlzhdrduBH83mbhM1klyKV+MVUPr8L9RDahRrigz912n8S9HWOh+U02ykVukdLc0q1/qKZd3r2rVqjIgIEBKKWX9+vVlSEhIqu2TPlWSHq8/NSOllC1btpR37tzJeLBSygcPHsgzZ87Ib775Rv76668pttu9e7dMSEiQCQkJsk+fPvKPP/6QUkr5888/y6+++kpKKeXDhw9lsWLFZHR0tAwNDZW2trYyNDRUhoWFSVtbWxkWFpapGHOKuPBweX/aNOlj7yCv1XSRIYv+kvEvXrxVn+fvhsvBy09L6wm7ZI0f98t5h3zl4+cxOopYySxSecpFnaHnEyNHjsTPz4927drxyy+/YGRkRIkSJYiIiMDW1lZ7Vvv06dNXlrds2UKdOnWoUqUKnp6egOZMvHHjxri4uODi4sKJEycAmDhxIp6enjg7OzNnzhwAOnXqxMaNmRs4bGFhQe3atdOcqKJ9+/YIIRBCUKdOHQIDAwHNWWxERITm8kRkJMWLF8fAwID9+/fTunVrihcvTrFixWjdujX79u3LVIw5hX7RopT6+msq7NxBwQb1CZkzh1vt2/Nk1+5MXV8HcC5XlBVD6rB9dENcrYsz+6AvjWZoBiiFq5GnOZJK6PnEokWLKFu2LO7u7lhYWODi4gKAqakpzZo1Y/fu3QBs3LiR7t27a5NoXFwcZ86cYe7cufz444+AJtEePHiQc+fOsWnTJsaNGwfAjBkzaNy4Md7e3owfPx4AV1dX7S8Cd3f3ZAchNWjQQCefMTY2ljVr1tC2raYc/5gxY7h69Sply5bF0dGRefPmoaenl2fKBiTHyNaWcgsWUH7VKvSLFuXeF1/g36cPz8+fz3SfNcoVZekgV3aPa0TjKiW0A5Sm771KSES0DqNX3pa625EPBQcHk7R88bBhw5g5cyZdu3ZlxYoVLFmyRLvt5ejNWrVqacsCxMbGMmbMGLy9vdHX18fX1zfFY1lYWGirMzZv3lz73HlWGDVqFE2aNKFx48YA7N+/H2dnZ44cOcKtW7do3bq1dlteV6huHWy3bOHJ9h2EzJnDnb79KNK+HSU/+5wCVpkbUWtf1ow/+tfC90EEC91vsuSYH6tO+NO3Tnk+alKR0mbGaXeiZCl1hp4PmZiYvDKUv2HDhvj7++Ph4UF8fDwODg7abS9LA7wsCwCa+i2lSpXiwoULeHl5EROT8p/fUVFR2sFHaZ2hL1y4ULsuIyV6AX788UdCQkKYPXu2dt2KFSvo3r07QggqVaqEra0t165dy7NlA14n9PUp2r0bFffvo8To0UQcccevfXse/vYb8ZGZr5tepZQp8/rU5NBnTenoVJbVJ+/QZKY7X/99iTuhz3T4CZSMUgk9H6pevTo3b958Zd3AgQPp168fQ4YMSXP/J0+eUKZMGfT09FizZo22DrmpqSkRERGvtPX19dX+gnh5hv766+U1+NGjR2vXlS1bNt2fZ+nSpezfv58NGzagp/ffP+ny5ctz+PBhAB48eMD169epUKECbdq04cCBA4SHhxMeHs6BAwdo06ZNuo+X2+gVLEjJsWOouG8vRdq3J3TJUm6914bwjZuQcZl/eqVCycLM6lUDjy+a0bt2ObaeC6T5LA8+2Xie6/cj0u5A0b2U7pZm9Us95fLuWVtby5CQEPns2TNpZ2cnExIStNuCg4OlsbGxDA8P165L+pRLSEiItLa2llJK6evrKx0dHaWTk5P86quvZKFChaSUUsbExMjmzZtLJycnOXv2bCmlpj7Mjh07MhVvcHCwtLS0lKamptLMzExaWlrKJ0+eSCmlbNeunQwKCpJSSqmvry8rVKgga9SoIWvUqCF//PFHKaWUQUFBsnXr1tLBwUHa29vLNWvWaPtetmyZrFixoqxYsaJcvnx5puLLrZ5fuiz9+38gfapWk7c6dpQRxzx10u+DJy/ktN0+0u67vdJ6wi45bNVZef5uuE76Vv5DKk+5qIFF+dQnn3xCp06daNWqFQBubm5s376dNWvW6OwY0dHRNG3alOPHj2OQyeHpStaQUhJx6BAPf51F7N27FGrcmFJffYmRDkYRP34ew8oT/qz4x58nL2JpUNGcUc0q0bCSOUIIHUSfv6U2sEgl9HzqwYMHnD59ms6dOzN27Fj27t3Lnj17qFIl8yMNX3fjxg2CgoJo1qyZzvpUdEvGxBC2fj2P/viThMhIir7fi5Jjx2Jgbv7WfUdGx7H+9B2Wet7mYUQ0jpZmfNysIm3sS6OvpxJ7ZqmErihKquLCw3m08A/CN2xAz8SEEiM/otiAAeglqZefWdFx8fzvXBB/HfPj9qNnVChRiI+aVqBrTUuMDPJWhcx3QSV0RVHSJdrPj4e/ziLS3R1DS0ssvvgc07ZtdXKpJD5Bsv/Kff7wuMnloKeUKmLEh41s6VfXWtWLyQCV0BVFyZBnJ0/yYMYvRF+/jknNmpSaOAGTGjV00reUkn9uhvLn0Zv8czMUU2MDPqhnzZAGNlgUUc+yp0UldEVRMkzGx/Pkf//j4bx5xIc8okjHjlh8Nh7DDDxSmpYLAY9ZfMyPvZeDMdDTo1tNS4Y3qUAli8I6O0ZeoxK6oiiZFh/5jNBlSwlbrqlWWXzwYMyHD0e/cCGdHcP/0TOWHvdji1cg0XEJtKpeipFNK+BqU1xnx8grVEJXFOWtxQYH83DOHJ7u2Il+iRKU/GQcRbt3R+hw6r/QyGhWnbzD6pP+PH4ei0v5ooxoUoHWdurJmJdUQlcURWdeXLzIgxm/8OLcOYyqVKHAR19i1aGRTo/xPCaOLV6BLPH0IzD8BdbmBfmwkS09a1lRsED+voGqErqiKDolpSRi/wFu/L6WU+WGUKJoPK69a2BboyR6OjyTjotP4IDPAxYf88M74DFFCxryQV1rBjawxsI0f95AfeuELoRoC8wD9IGlUsoZr20fDPwKvKxBukBKuTS1PlVCV5TcR0rJkydP8PX1pVq1ahgJQ87/sRffJ6WICIumSEkTarQoR/UGZTA00t2lGCkl/94JZ4mnHwd8HmCop0cX57IMa1yBqqVNdXac3OCtEroQQh/wBVoDgcBZoK+U0idJm8GAq5RyTHqDUgldUXKPbt26cenSJe7du0dCQgLR0dFMmzaNr7/+GoCEBInf+RC8D93lwe2nGBU0wKGJJY7NrShk9vaDk5Lyf/SM5f/cZrNXAFGxCTSuXIKhjWxpWlm3fx3kVG+b0OsDP0gp2yQufw0gpZyepM1gVEJXlDzrypUrBAQE8NNPP3Hv3j0KFSrEhQsXXqlu+VLwrSd4H7yL34UQ9PQFVeqUxrllOcwtdfsoYvizGNafucvqk/48eBpNxZKFGNrIlu41rTApkHdHoL7tJNGWQECS5cDEda/rIYS4KIRwE0KUS2a7oii5UHR0NFu3buWDDz6gY8eOPHv2jJUrVyabzAHKVDSj3UhH+v9YD/uGZbl59gEbp55h53xvAq6GZXpKvNcVK1SA0c0r4flVC+b2dqZgAQMm/e8yDWYc5tf913jwNCrtTvKY9Jyh9wTaSimHJS4PAOomPRsXQpgDkVLKaCHER0BvKWWLZPoaAYwAKF++fK07d+7o7pMoiqJznp6ejBgxgipVqrBgwQKmT5+OEIKFCxemu4+oyFguHwviokcgL57GYG5ZGOfW5ajsWgp9A91NySCl5Kx/OMuOa66zG+gJ2juWYUhDW5zLFdXZcbJbll9yea29PhAmpTRLrV91yUVRcq7Hjx8zYcIEdu/ezfz58+nWrRtCCP744w/69u1LsWLFMtxnfGwC18/cx/tQAOHBzyhkVgCnFuWwb1wWo4KpTwSeUXdDn7PixG22eAUSGR1HzfJFGdLQlnYOpTHUz93z+rxtQjdAc1O0JZqnWM4C/aSUV5K0KSOlDE583w2YIKWsl1q/KqErSs4jpcTNzY1PP/2ULl26MH36dMzMUj03y9Qx7vqE4X3wLoHXwjEw0seuYRlqtChHkRImOj1WRFQsW/8NZOUJf/xDn1OqiBED6lnTt055zAvr9mbtu6KLxxbbA3PRPLa4XEr5sxBiCpqZM3YIIaYDnYE4IAz4WEp5LbU+VUJXlJzl7t27jB49Gj8/PxYvXkzDhg2z/JghARF4H7rLzbMPkVJS0cUC51blKWVbRKfHSUiQePg+ZMU//njeeEQBAz261CjLoAY2OFjq9hdWVlMDixRFSVF8fDy///47P/30E+PHj+fLL7+kQIEC7zSGyPAoLroHcsXzHjEv4ihTyQznVuWxcSqh80cRbzyIYOUJf/4+F8SL2HhcrYsxqIENbXPJ5RiV0BVFSdb58+cZMWIEhQsX5q+//tLpjFWZERMVx9V/grlwOICIsCjMLExwblmOqvXLYKjjRxGfvIhli1cAa07d4U7ocyxMjehf15q+dcvl6FGoKqErivIGKSWNGjXiww8/ZMiQITlqvs+E+ARunQ/B++BdHt6JwLiQIQ5NLXFsZkXBIrr96yEhQXLUN4SVJ/w56huCob7m6ZiB9W1wKV80R30voBK6oii5lJSS4JuPOX8wAP9Lj9DX16Nq3VLUaFWe4mV0V773Jb+QSNacuoObVyAR0XHYlSnCwPrWdHYum2OKgqmErihKrvf4wXO8Dwdw7WQw8bEJWDuY49y6PJZVdH8W/Sw6jm3eQaw5eYdr9yMwNTagZy0rBtSzpkLJ7J18QyV0RVHyjBcRMVw+FsQlj0BeRMRSolxhnFuVp5KrBfo6vqkppcTrTjhrTt5h7+VgYuMljSuX4IN61rSsZoFBNtxEfduh/4qiAB4eHnTs2DFD+6xcuZJ79+6lu72UknHjxlGpUiWcnJw4d+5csu2aNWtG1apVcXZ2xtnZmYcPHwJw7NgxXFxcMDAwwM3NTdve29ub+vXrY29vj5OTE5s2bcrQ58hJTEwLULuDLQN/bkCz/lWJj03g0Aof1n57kvMH7hL9Ik5nxxJCUNumOPP71uSfiS344r0q3HwYyUdr/qXRL+7MPeTL/Sc5p8RAzrgopCh51MqVK3FwcKBsOufh3Lt3Lzdu3ODGjRucPn2ajz/+mNOnTyfbdt26dbi6vnqiVr58eVauXMmsWbNeWV+wYEFWr15N5cqVuXfvHrVq1aJNmzYULVo0U58rJzAooI99Y0vsGpblzuVQvA/d5cTfNzm75zZ2jcpSo0U5TIvr7mkVC1NjxrSozMimFTl87SHrTt9l7qEb/H7kJi2rWfBBPWsaVdL9Y5YZoRK6oiRj8uTJFC9enE8//RSASZMmYWFhQWRkJD179uTy5cvUqlWLtWvXIoRgypQp7Ny5kxcvXtCgQQP++usvtm7dipeXF/3798fExISTJ09iYpL6SMjt27czcOBAhBDUq1ePx48fExwcTJkyZdIVt42NDcAbhbOSPo5YtmxZLCwsCAkJydUJ/SWhJ7BxKoGNUwke3nmK96EALh4J5OKRQCq5lMS5dXksrHU3UMlAX4829qVpY1+au6HPWX/mLlu8Ajjg84DyxQvSr255etWyypaRqOqSi6IkY+jQoaxevRqAhIQENm7ciJWVFefPn2fu3Ln4+Pjg5+fHP//8A8CYMWM4e/Ysly9f5sWLF+zatYuePXvi6urKunXr8Pb2xsTEhPHjx2svkyR9zZihmTMmKCiIcuX+K1ZqZWVFUFDQmwECQ4YMwdnZmalTp2aoguGZM2eIiYmhYsWKmf16ciwL6yK896E9A36qT40WVvhfDmXLdC/+99s5bl98hEzQ7T3D8uYFmdiuGie+bsG8Ps6UNjNmxt5r1Jt+mDHrz3Hi1iOdVZdMD3WGrijJsLGxwdzcnPPnz/PgwQNq1qyJubk5derUwcrKCgBnZ2f8/f1p1KgR7u7uzJw5k+fPnxMWFoa9vT2dOnV6o985c+boJL5169ZhaWlJREQEPXr0YM2aNQwcODDN/YKDgxkwYACrVq1KsfxtXmBa3JiGPSvj2sEWn+P3uHgkgD1/XKRoqYI4typH1bqlMdDhQCUjA326OFvSxdmSGw8iWH/mLlv/DWTXxWBsSxSiT+1y9KhlRYksPmtXCV1RUjBs2DBWrlzJ/fv3GTp0KABGRv/9QOrr6xMXF0dUVBSjRo3Cy8uLcuXK8cMPPxAVlfyNsvHjx+Pu7v7G+j59+jBx4kQsLS0JCPhv+oHAwEAsLd+cfuDlOlNTU/r168eZM2fSTOhPnz6lQ4cO/Pzzz9Srl2rtvDzDyMSAmq3L49TCilv/PsT7UAAe665zeocfDk2tcGxqiYmpbgcqVS5lyved7JnQthp7LgWz4cxdpu+9xqwD12ljX5p+dcpTr4J5llxrVwldUVLQrVs3Jk+eTGxsLOvXr8fT0zPZdi+Td4kSJYiMjMTNzY2ePXsCmoQbERGhbZvWGXrnzp1ZsGABffr04fTp05iZmb1x/TwuLo7Hjx9TokQJYmNj2bVrF61atUq135iYGLp168bAgQO1seUn+vp6VKlTmsq1S3HP9zHnD93l7K7bnNt/h6r1NDMqFSut24FKxob6dHexoruLFTceRLDhTABbz2nO2r9uV42Pmur+kpdK6IqSggIFCtC8eXOKFi2Kvn7Kf54XLVqU4cOH4+DgQOnSpaldu7Z22+DBgxk5cmS6b4q2b9+ePXv2UKlSJQoWLMiKFSu025ydnfH29iY6Opo2bdoQGxtLfHw8rVq1Yvjw4QCcPXuWbt26ER4ezs6dO/n++++5cuUKmzdv5tixY4SGhrJy5UpA8wSOs7Nz5r+gXEgIgWXVYlhWLUZY8DMuHA7g+sn7+Hjew8apBM6tylG2su4HKlUuZcrkTnZ81bYqey8HU9fWXKf9v6QGFilKChISEnBxcWHLli1Urlw5u8NRssjzpzFcOhrIZY8gop7FYmFtinOr8lR0KYleDqy+qAYWKUoG+fj4UKlSJVq2bKmSeR5XsEgB6naqwMDpDWjaryrRL+I4sOwKa787hfehu8RE6W6gUlZTZ+iKoihJyATJ7YuP8D50l+CbTyhgrI9dY0ucmlvpdKBSZqkzdEXJhPnz51O9enX69+/Ptm3bmDJlSqrtmzVrRkZOUvz9/Vm/fr12+dKlSwwePDiz4XLt2jXq16+PkZHRGyNFk7p9+zZ169alUqVK9O7dm5iYGABmz56NnZ0dTk5OtGzZktcncX/69ClWVlaMGTMmuW7zDKEnqOBcku5f1KLnBFfK25tz4dBd1n57koPLrxByNyLtTrKJSuiKkoI//viDgwcPsm7dOmbOnMmoUaN02v/rCd3R0ZHAwEDu3r2bqf6KFy/O/Pnz+eKLL1JtN2HCBMaPH8/NmzcpVqwYy5YtA6BmzZp4eXlx8eJFevbsyVdfffXKft999x1NmjTJVGy5VSnbIrQZ7sAHU+vj2MyK2xcesXnaWbbNOY//Jd0PVHpbKqErSjJGjhyJn58f7dq145dffsHIyIgSJUoQERGBra0tsbGxgOasNenyli1bqFOnDlWqVNE+5ujv70/jxo1xcXHBxcWFEydOADBx4kQ8PT1xdnbWPs7YqVMnNm7cmKmYLSwsqF27NoaGhim2kVJy5MgR7aOLgwYNYtu2bQA0b96cggULAlCvXj0CAwO1+/377788ePCA9957L1Ox5XZFSpjQ6P3KDJregPrdKvL4wXN2L7zIhqln8PnnHnGx8dkdIqASuqIka9GiRZQtWxZ3d3csLCxwcXEBNM+VN2vWjN27dwOwceNGunfvrk2icXFxnDlzhrlz5/Ljjz8CmkR78OBBzp07x6ZNmxg3bhwAM2bMoHHjxnh7ezN+/HgAXF1dtb8I3N3dky0T0KBBg0x/rtDQUIoWLYqBgeaJ5ZRKCyxbtox27doBmqd9Pv/881Qv4+QXRgUNcWljzYCf6tNqcHX09AXua66x+psTeO25TVRkbLbGp55DV5Q0BAcHU7JkSe3ysGHDmDlzJl27dmXFihUsWbJEu6179+4A1KpVC39/fwBiY2MZM2YM3t7e6Ovr4+vrm+KxLCwstOV2mzdvjre3t+4/UBrWrl2Ll5cXR48eBTSXntq3b68teaCAvoEeVeuVoUrd0gReD8f74F1O77jNv3vvUK1+GWq0LEfRUgXfeVwqoStKGkxMTHjy5Il2uWHDhvj7++Ph4UF8fDwODg7abS9LA7wsCwCa0aGlSpXiwoULJCQkYGyc8pMSUVFR2sFH7u7u2jP3pAoWLMiJEydYuHCh9pfJnj170lWi19zcnMePHxMXF4eBgcEbpQUOHTrEzz//zNGjR7Wf5eTJk3h6evLHH38QGRlJTEwMhQsX1hYUy8+EEJSrVpxy1YoTGhSJ9+EAfE7c47JnELZOJXBuXZ4yFc3e2bykKqErShqqV6/O2rVrX1k3cOBA+vXrx3fffZfm/k+ePMHKygo9PT1WrVpFfLzmeuvrZQEAfH19tb8g0jpDHz16NKNHj87QZxFC0Lx5c9zc3OjTpw+rVq2iS5cuAJw/f56PPvqIffv2YWFhod1n3bp12vcrV67Ey8tLJfNkmFsWpuXA6tTrUoFLHoFcPhrE7QuPKGVbBOdW5angXCLLByqpa+iKkoYmTZpw/vz5V8qg9u/fn/DwcPr27Zvm/qNGjWLVqlXUqFGDa9euUaiQpmaIk5MT+vr61KhRQ3tT1N3dnQ4dOmQqzvv372NlZcXs2bP56aefsLKy4unTp4CmpMDLSzm//PILs2fPplKlSoSGhvLhhx8C8OWXXxIZGUmvXr1wdnamc+fOmYojvytkZkS9LhUZNL0hTfpU4UVkLPuXXGbt5FNcOBKQpQOV0jWwSAjRFpgH6ANLpZTJ/noWQvQA3IDaUspUH8hVA4uU3OSTTz6hU6dO2iJYbm5ubN++nTVr1ujsGNHR0TRt2pTjx49rb1oquV9CguT2hRC8DwZw3+8JRgUNaNK3ClVql85Uf6kNLErzX40QQh9YCLQGAoGzQogdUkqf19qZAp8Ayc+XpSi52DfffKOdCm7s2LHs3buXPXv26PQYd+/eZcaMGSqZ5zF6eoKKNS2oWNOC+35P8D54lyIlUi/SlllpnqELIeoDP0gp2yQufw0gpZz+Wru5wEHgS+ALdYauKIqie2879N8SCEiyHJi4LukBXIByUsrdmY5SURRFeStvfVNUCKEHzAY+T0fbEUIILyGEV0hIyNseWlEURUkiPQk9CCiXZNkqcd1LpoAD4CGE8AfqATuEEG/8SSClXCyldJVSuiYdqKEoiqK8vfQk9LNAZSGErRCiANAH2PFyo5TyiZSyhJTSRkppA5wCOqd1DV1RFEXRrTQTupQyDhgD7AeuApullFeEEFOEEOpBVUVRlBwiXc9HSSn3AHteWzc5hbbN3j4sRVEUJaPUSFFFUZQ8QiV0RVGUPCLb5hQVQoQAd9JsmLwSwCMdhqMrKq6MUXFlXE6NTcWVMW8Tl7WUMtnHBLMtob8NIYRXSiOlspOKK2NUXBmXU2NTcWVMVsWlLrkoiqLkESqhK4qi5BG5NaEvzu4AUqDiyhgVV8bl1NhUXBmTJXHlymvoiqIoypty6xm6oiiK8pocl9CFEG2FENeFEDeFEBOT2W4khNiUuP20EMImybavE9dfF0K0yQlxCSFshBAvhBDeia9F7ziuJkKIc0KIOCFEz9e2DRJC3Eh8DcpBccUn+b52vL5vFsf1mRDCRwhxUQhxWAhhnWRbdn5fqcWVnd/XSCHEpcRjHxdC2CXZlp0/j8nGld0/j0na9RBCyKRFDHXyfUkpc8wLzRR3t4AKQAHgAmD3WptRwKLE932ATYnv7RLbGwG2if3o54C4bIDL2fh92QBOwGqgZ5L1xQG/xP8WS3xfLLvjStwWmY3fV3OgYOL7j5P8f8zu7yvZuHLA91UkyfvOwL7E99n985hSXNn685jYzhQ4hqaQoasuv6+cdoZeB7gppfSTUsYAG4Eur7XpAqxKfO8GtBRCiMT1G6WU0VLK28DNxP6yO66slGZcUkp/KeVFIOG1fdsAB6WUYVLKcDSzTbXNAXFlpfTE5S6lfJ64eApNuWjI/u8rpbiyUnrieppksRDw8qZctv48phJXVkpPngCYCvwCRCVZp5PvK6cl9DRnR0raRmoqQT4BzNO5b3bEBWArhDgvhDgqhGiso5jSG1dW7JvVfRsLzUQop4QQXXUUU2bi+hDYm8l931VckM3flxBitBDiFjATGJeRfbMhLsjGn0eR8uxuOvm+1Gy0WS8YKC+lDBVC1AK2CSHsXzuDUF5lLaUMEkJUAI4IIS5JKW+9ywCEEB8ArkDTd3nctKQQV7Z+X1LKhcBCIUQ/4FtAp/cXMiuFuLLt51H8N7vb4Kw6Rk47Q09rdqRX2gghDAAzIDSd+77zuBL/hAoFkFL+i+baWJV3GFdW7JulfUspgxL/6wd4ADXfZVxCiFbAJDQTtURnZN9siCvbv68kNgJdM7nvO4krm38eU5vdTTffV1bcHHiLmwoGaG422fLfTQX719qM5tWbj5sT39vz6k0FP3R3E+Zt4ir5Mg40N0uCgOLvKq4kbVfy5k3R22hu8BVLfJ8T4ioGGCW+LwHcIJkbS1n4/7Emmh/yyq+tz9bvK5W4svv7qpzkfSfAK/F9dv88phRXjvh5TGzvwX83RXXyfb31h9D1C2gP+Cb+452UuG4KmrMSAGNgC5qbBmeACkn2nZS433WgXU6IC+gBXAG8gXNAp3ccV2001+OeoflL5kqSfYcmxnsTGJIT4gIaAJcS/3FfAj58x3EdAh4k/v/yBnbkkO8r2bhywPc1L8m/b3eSJLBs/nlMNq7s/nl8ra0HiQldV9+XGimqKIqSR+S0a+iKoihKJqmEriiKkkeohK4oipJHqISuKIqSR6iEriiKkkeohK4oipJHqISuKIqSR6iEriiKkkf8Hzi4fU2UiZHnAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Bonus cell. This is the cell Stephanie used to make a figure\n", "# in the PDF lecture notes. It illustrates the first 5 steps\n", "# taken to solve the test equation with a rate constant of 2\n", "# and step size of 0.1.\n", "params = np.array((2,))\n", "\n", "y0 = 1\n", "tru_ts = np.arange(0,0.41,0.01)\n", "tru_y1 = y0*np.exp( -params[0]*tru_ts)\n", "\n", "plt.plot( tru_ts, tru_y1 )\n", "\n", "# Illustrate forward Euler for a few steps\n", "# of size 0.1\n", "ts = np.arange(0,0.5,0.1)\n", "y = np.zeros((1,ts.size))\n", "y[0,0] = y0\n", "for i in range(1,ts.size):\n", " h = t[i]-t[i-1];\n", " yprime = degradation_model(t[i-1],y[0,i-1],params);\n", " y[0,i] = y[0,i-1] + h*yprime;\n", " plt.plot( (ts[i-1],ts[i]), (y[0,i-1], y[0,i]) )\n", " if i == 1:\n", " offsets = (0.01,0.1)\n", " else:\n", " offsets = (-0.1,-0.1)\n", " plt.annotate( \"yhat=%0.3f\\nf(yhat)=%0.3f\"% (y[0,i-1],yprime), (ts[i-1],y[0,i-1]), (ts[i-1]+offsets[0],y[0,i-1]+offsets[1]), \n", " xycoords='data', arrowprops={'arrowstyle':'->'})\n", "plt.legend(['True Solution (y)','Approximate Solution','Approximate Solution','Approximate Solution','Approximate Solution']);\n", " " ] } ], "metadata": { "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": 4 }