{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forward Euler, continued\n",
    "\n",
    "First, revisit the simulations with large time steps. On the board, we determined that, in order for the simulation to be constantly shrinking, that the time step should be < 1. \n",
    "\n",
    "Let's try that with our single-state degradation model. Here is copy-pasted code from Tuesday."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABERUlEQVR4nO29e5Qdd3Xn+9n9VrfULfXj6N2WbEnu0zZ+QFuEwIXBNmASbEOAYLFmMFwyXrNuCEmYOwGGdUNChlxyuTeQzGXNXAcIJjA2xsDYBgdjjAkkwaA2fqpPy3pYtt511C21pJb6/bt//KpOV7fO6T6Petfvs1avPqdOnapfV++q/Xvsvb+ilMJgMBgM6aUu7AYYDAaDIVyMIzAYDIaUYxyBwWAwpBzjCAwGgyHlGEdgMBgMKach7AZUQ3d3t9qyZUvYzTAYDIZY8dRTT51SSvUs3h5LR7BlyxYGBwfDbobBYDDEChF5udh2MzVkMBgMKcc4AoPBYEg5xhEYDAZDyjGOwGAwGFKOcQQGg8GQcjxxBCLyVRGxROSFEp+LiPytiOwXkedE5NWuz+4UkX32z51etMdgMBgM5ePViOBrwC1LfP52YLv9cxfw3wBEpBP4NPBaYCfwaRFZ41GbDAaDwVAGnjgCpdTPgNEldrkd+LrSPAmsFpH1wNuAx5RSo0qp08BjLO1QauLp79/Nr779eb8OHz+efwAuLPVvSxHjI/p6GDQv/wKOPxt2K6KBUvDrf4Cp8bBb4htBrRFsBA673h+xt5XafgkicpeIDIrIYD6fr64VuYfI7PkKE9Oz1X0/SZx+Gb7zYRj8StgtiQaDX9XX48zh5fdNAw99BH74ybBbEQ2OPqWvR4I7CrFZLFZK3a2UGlBKDfT0XJIhXRZdW69lszrBY88VTa5LF1Zu4e+0Y+2xf5vrwfRFGDkA1pDuDaedk8m3jaAcwVFgs+v9Jntbqe2+sOnK11Avin958hd+nSI+WEP27+Qad0UUHONQuO2IAvm9gIKLp+H8ybBbEz4psI2gHMFDwAfs6KHfAMaUUseBR4G3isgae5H4rfY2X6hb2w/AxNEXeOlUcuf7ysIx7lMvwsxUuG0Jm5lJGNmvXxvHuPAaJPjhVzYp6DR5FT56L/AL4EoROSIiHxaR/yAi/8He5RHgILAf+DvgfwNQSo0CfwHstn8+Y2/zh64rUHWN9NUf4b7dr/h2mlhg5UDqYW4GRg+E3ZpwGdmvr4PUmwcf6Gsg9fbr5D78ysa5V8YtHVSQQLyKGtqllFqvlGpUSm1SSn1FKfXflVL/3f5cKaV+Xyl1hVLqVUqpQdd3v6qU2mb//L0X7SlJfSPSvYPfXHWS7zx1hKmZOV9PF1lmZ/RI4PI36fdpv9mdv//yN+nrMpfyYIL8MGT6oa3H2Mb4iHYAzr2ST+b1iM1isWdk+tguRzl1forHcymd/zz9EsxOQvY2uxecTOMuGysHdQ36esxMwOlDYbcoXKwcZLL6J+224Tz4r363/p3Q65FCR5BlxfgRtrUr/sevUjo95Ex/bLgOOi830yFWDjqvgPXX2O9TfD0mzsLYYcj0QU9Wjw7mUjpyBtdo8c3Q0pFY20ihI9ALxh/um+Kf95/i8OiFkBsUAlYOEOi+0vT6QN/cmSz09NnvU3w98sP6d6ZfX5Op89oxpBVrCJo7oH2DviYJtY0UOoIsALdkTiPA/YMpNHJrCDq3QlOrNu7Rgzp2PI1MXdBTQZl+aGqDNVsS2+srC+dvz2QLnaakPvzKwpkmE7E7TcnMrUifI1i9BRpWsOb8ft60o4f7Bw8zM5uyoa81PH+TZ7KA0oukaeSUHTNvdxB0r2841CaFijUMjW3Q0aunhyCxC6TLotS8IwBtGxNjcO5EuO3ygfQ5gro66LkSrBx37Ozl5NlJnthbZcmKOOLEzLuNG9Lb63P+brdjHNmX3twKa0g7gLo6PSfevim9tnHuBEycWdRpIpEjxvQ5AijM9d3YlyGzqpn70rRofGofqNn5+fDOy6G+KZHGXRbWENQ366ky0AukczPzCWZpw8rpa+CQ6Uu3bcD8yMi5Lgl0jCl1BFk4f4LGyTO8d2ATT+y1OD6WkjnyxT3g+gbo3pFI4y4LKwc9O6DOTqBKcK9vWcZP6Zj5jNsRZCH/os49SRuL75W2LmjLJPJeSakjmJ8Oed9AL3MK7t99JNw2BYU1pGPmu7bNb0tz5JCVm7cHgO7t6c2tKDz43I6gX+ecnH4pnDaFiZXTSXVt3fPbnAXjhJFSRzDf6+vtauUN27q5f/Aws3PJiwa4hPwwdG2Hhqb5bZmsDhGcOBteu8JgYgzOHl344Gto1k4yn8IFY3foqEMmudMhy5LPLbQN0NcmvzdxuRXpdATtG6C5vWD4u3b2cvTMRX6+LwWLxk7MvBvnxs/vDb49YWIVefBBYnt9y2INQctqWLVuflv3lYCkzxHMzS2MrnPIZGF6HMaSta6YTkdQiAnWxv2W/rV0tTVxb9IXjafG7Zj5RY6gkEiVsoef8/c6f79DJgujL+kcgzThjpl3aGpNZ27F2Cv6gV/MNiBxjjGdjgAWJIc0NdTx7tds4vGchXVuIuyW+Udh6L/IEay+DBpbE2fcy2LloGkldGxeuL2QW5GiEZJSxUeLkOiM2pIsXih2SGinKcWOoH+B8MYdN2xmZk7xwFMJXjQuNRVSV6cNPG2JQ9aQ/rvrFt0GhWCCFK0TnDuu10wW2wbYI6QDOgclLRQcwaIRQUu77jgkzDGm1xEs8uyX96zktVs7+dbuw8wlddHYGoKGFj3UX0wae3354eI94DVbdW5Bwnp9S1IsYsghk8LcCiunk+laOi79LJNNXCchvY6gSK9v185eXh65wC8OJlN8QsfMXzkfM+8mk9Wjo4QKb1zC+TyM54v3gOsbdG5Bmhyj87f2lJgacu+TBqwiEUMOmayeNkxQboVXCmW3iMheEdkvIp8o8vkXROQZ++dFETnj+mzW9dlDXrSnLFb2QGv3gl7fLVevo2NFY3IXjRdnjbpJW12ZfImhv0NPynIrrJxOlmrruvSzrm069yQtI6TZGf2gX8o2Zqd0scaEULMjEJF64EvA24F+YJeILOhmKaX+WCl1nVLqOuC/At91fXzR+UwpdVut7amIRYlULY31/M6rN/KjPScZOZ+w+dCLZ+DcsSV6OSnr9ZVaDHTIZOHsET1vngZKLRSDzjnp2pYe2xg9qB/0S9kGJMoxejEi2AnsV0odVEpNAfcBty+x/y7gXg/OWzuZ/kuEN3bt7GVqdo7v/vpoiA3zgWLJQm5WrbeFN1Jys1tDsGINrFxb/PM05VbMzc3LU5YiTdnn+SXWS0BPryYst8ILR7ARcBf1P2JvuwQRuQzYCvzEtblFRAZF5EkReWepk4jIXfZ+g/m8R4lfmb5LhDd2rF3Fq3tXc+/uV1BJqjvurjNfDJF0LRg7yULumHk3Cez1lWTsFZi+UNo2QF+r04d0LkrScQs3FaNxhS7WmKBp1KAXi+8AHlBKudXBL1NKDQDvB74oIlcU+6JS6m6l1IBSaqCnp8eb1hR6fQsjAHbt7OVgfpzdh057c54oYOWgaRV0bCq9T4KFNxawuM58MTo26xyDNDjG5abJYD63Ig0jJLdwUykSNkLywhEcBdwZOZvsbcW4g0XTQkqpo/bvg8BPges9aFN5lEgO+e1r1rOquSFZi8ZWTo+ASvWAQS+CTZxJpPDGAs4eg8mxS7NG3RR0K1IwIihkWJfoAUOiSzBfwlJBFQ6ZLIwcgOlkJKB64Qh2A9tFZKuINKEf9pdE/4hIH7AG+IVr2xoRabZfdwOvB4K781ashvaNlxh3a1MDt1+/gUeeP87YhenAmuMrSy0GOqRlOqScHjAkrtdXEiunR0At7aX36UxJbsXMpH7Al3OvqFktYpQAanYESqkZ4CPAo0AOuF8ptUdEPiMi7iigO4D71MKJ9ywwKCLPAk8An1NKBWtpJQqM7drZy+TMHN97OgGZxufzcGGkvAcfJL/y5nLrJQ6Zfp1rMH7K/zaFyXLTZKBzT3quTL5tOMJN5dgGJKaj0ODFQZRSjwCPLNr2p4ve/1mR7/0r8Cov2lA1PX3w0s9hbnZBotVVGzq4ZlMH9+0+zJ2/uQVZakol6pT74GvrtoU3Et7ryw/DynXQ2rn0fu4CY1v/F//bFQazM1qvettNy++b6YdDP/e/TWFS7mix8wqoa0yMI0hvZrGDI7wxeqnwxh039DJ84hxPHz4TfLu8pFzjhnRMh5QzTQaJ6/UVZbmYeTeZrNZvuHjG92aFRjHhpmIkLLfCOIIl5sVvu24DrU318dc0toZgRadWW1oOp45KwoQ3ChTqzJfhCFau1fX5kzxCKlWKuxhpmDq0cvoB7xZuKkWCdCuMI1giOWRlcwO3XbuBh589zrmJGC8aO3KM5UxvJVR4o8CZQzBzsTxHkIbcCidmfqmIIYc0BBOUO1oEbRtnXobJ8/62KQCMI2hqW1J4446dvVycnuXBZ44F2y6vUKp0lc1iJL0EcyXTZGCLt+eSm1thDenkqMYVy+9byK1IqG1MjesHeyW2AYnIrTCOAJacF792Uwd961Zx3+6Y9pDPHoXJs+U7goQKbxQoVNksowcM+rpNjOl6/Umkkk5CQdkvobZRSripFAVHEP8Ro3EEsKTwhojw/tf28sLRszx/JIYFyJaqM1+MlnZdhz2p0yFWDjp6oXlVefsneTpkeqK8mHk3PX3Jtg0of0SwZovW90jA9TCOAPQ/fgnhjduv20hLYx33xnFUUMlioEOSI4fKiZl3k+SM2pEyY+bdZPrhwimdm5I0rFxp4aZiOLkVCegkGEcAywpSd6xo5LdetZ6HnjnG+GTMxCisnK4sulzMvJsECm8AMDutY+YrefC1denooSQ6gkp7wJDsEZI1BN07igs3lSIhwQTGEQB0bV9WeOP9O3s5PznDD56L2VxxpT1gsHMrpuD0pbkVsWbkAMxNV/bgg+SOkKwhnRTVWbTOY3FKFGpMBE5F2krIZPX60cV4F6g0jgB0zHDnFUve7K+5bA3bMivjNT00N6sjGqoxbkher2+5OvOlKKJbkQisYejeXl7MvMPKjM5JSZptXDy9tHBTKRISZWccgcMyvT4R4Y4bNvP0K2cYPnE2wIbVwOlDOma+kvUB0MPjhAlvAPrvkTr776uAnj5dr//My/60KyysocptoxA5lDTbWEa4qRQJibIzjsChDOGNd796E031ddz3q8Ml94kU1cwBg67D3rk19sZ9CYWY+ZbKvpfEUhOT5yuLmXfjOIIk5VYU6nFV6Bg7Nmmdj5jbhnEEDmUIb6xpa+KWq9fx3V8fYWJ6tuR+kaHSmHk3CVkEW0A16yUwf/2S5BgdO6/memSyOjflbILkXK2cTpbr2Lz8vm4SMkIyjsChzF7fHTs3c3Zihkeej8GicT4Hqy+D5pWVf9cR3iiSWxFLpid0gbVqesAt7Tr3IEkLpOVWpC1GQubFF+Ak1lVTZTgByn7GETiUKbzxusu72NLVGo/pIafGUDU4whunkiG8wakXQc1V9+CDRPT6FpAfhoYV5cfMu0nIvHgBpeDknhpsox8ujmrtiphiHIFDXT307Fi21yci3LGzl18dGmW/FeFiUzNTdsx8hXOeDklLpCpMk1V7s/fp6zkb4+KDbqwhbe+VxMw7tHZqPYek2MZ4Xj/Iq+40xd8xeuIIROQWEdkrIvtF5BNFPv+giORF5Bn75/dcn90pIvvsnzu9aE/VlDkv/u5Xb6KhTvhWlENJRw/obOlqjbtr27K5FbHCiZnvqiBm3o2TWzF60Nt2hUUto0VIVs2harLv3SQgmKBmRyAi9cCXgLcD/cAuESlmYd9SSl1n/3zZ/m4n8GngtcBO4NMisqbWNlVNmcIbPauaeUv/Wh546giTMxFdNK5lDhhs4Y3tsTbuBVg5HTZa31jd95OUW3FhVCdBVWsbYOdW7E1GbkW10XUObT3Q2hVr2/BiRLAT2K+UOqiUmgLuA24v87tvAx5TSo0qpU4DjwG3eNCm6qgga3LXzl5OX5jmR3tO+tyoKrGGQer1w7xanBLMSaDaiCGH7h06ByEJC6T5KmPm3WSyOkflzCFPmhQqVk4nya3MVPf9gm5FfG3DC0ewEXCvnB6xty3m3SLynIg8ICJOjFa530VE7hKRQREZzOd9WpSpYBHsDdu62bRmRXTLU1tDehqk0ph5N2XkVsSCyXNaaKcWR9C4QucgxLjXV6DSirTFSMB0SIFKhJtKEfPciqAWix8GtiilrkH3+u+p9ABKqbuVUgNKqYGenjIkF6uhAuGNujrhfQOb+Zf9I7w8EsEHpZWrfs7TwVkEi3vYZC0x826SUoLZykFzO7QX7XOVR1JyK5SqfbQI2jamzsHYEW/aFTBeOIKjgDsLY5O9rYBSakQp5QSkfxl4TbnfDZS6OvtmL8+43zuwmfo64b7dEQslnb5Yfcy8m6T0+mpdL3HI9OtF+OmJ2tsUJk4noZYecPNKWN0bf9sYO6If4NVG1znE/F7xwhHsBraLyFYRaQLuAB5y7yAi611vbwOcq/Uo8FYRWWMvEr/V3hYeFcSLr+to4c1XZvj24BGmZyO0aJbfC6jaH3xJEd6wcjpmfvWW2o6TyepchFMvetKsUFCqMl3epUhC9nmtC8UOMQ8hrdkRKKVmgI+gH+A54H6l1B4R+YyI3Gbv9lER2SMizwIfBT5of3cU+Au0M9kNfMbeFh4VCm/s2rmZU+cneTwXoUVjLxYDwSW8EfebfUjfqHU1mnsSSjCft2qLmXeTyeqEwzjnVjjBELVOo65YA6s2xNY2Grw4iFLqEeCRRdv+1PX6k8AnS3z3q8BXvWiHJ7g9+8o3Lbv7m3b0sK69hXt/dZhbrl6/7P6BYA1BfZNe3KyVTD8c/KfajxMm1jBsu6n243RdoXMRYtrrA6ovxV2MTL/Wdxg5UPvUSlhUI9xUihjnVpjM4sVU2OtrqK/jd2/YzM/25Tly+oKPDauAQsy8B36+p0/XaY+r8MaFUTh/ovYeH+gchO6Y51Z4ETHkkIRSE15Nk4Edbr1X64DEDOMIFrNyrR7mVWDcvzuwCYD7ByMSMeBFFIRD3AuMeTUH7BDjXh+g297apZOgaqWQWxFTx+gIN1VbdmQxmSzMTOiQ65hhHMFiCskh5Rv3pjWtvHF7D/fvPsxM2IvGE2dh7LC3vRyI78PPq4ghh0wWzryi6/nHES9i5h0aW2xlv5jaxulD+sFt7hXjCIpSRXLIrp29nDg7wT+9GHIFwkLMvEc9YEd4I6aLYDpmvgPaN3hzvMLUYWndisiilK3L69GDD+zpkBjbBnh3rxSmyuJ3PYwjKEZPX8XCGzdlM3SvbObesMtTe90DjrvwRi115osR417ffMy8l46gX+esTF/07phBUYtwUzGa2rT+RwxtwziCYlQxL95YX8d7Bzbxk+GTnBgLMeHIykFjqxZS8YpMn67XHrf0+ULMvIcRLau36JyEODrGWktxFyPTF9/cCmuoeuGmUsQ0t8I4gmJU2eu744bNzCn49mCIowJHkLzWmHk3cRXeOH9SRzt5NfQHO/v8ylj2+qrW5V2KOGfUehlU4ZDJwsg+rQcSI4wjKEaVwhuXdbXx+m1d3Lf7MHNzIfWea60zX4y4Tod4PU3mENNen46Z36Cj4ryi83KdsxI325iZ0g9sP2xjbgZG9nt7XJ8xjqAUVYYJ3nFDL0fPXOTn+0/50KhlGB+Bccsf44b4LYJ5vRjokMnq3IQL4SbBV4yXMfMO9Y06jDRutlGrcFMpnOsbs/LtxhGUopAcUlk46FuvWsua1kbu+1UI5am9zBp1E1fhDSun297W7e1x41hqYm7Wli712DYgnsEEfo0Wu7drHZCYXQ/jCEpRpfBGc0M973nNJh4bOkn+3OTyX/ASL7NG3YjoBcaYGbcnpbiLEccCY17HzLvp6dN6DxNnvT+2X1i52oWbitHQrEuRxOxeMY6gFDUsgr3vhl5m5hTf+XXAmcbWELR06NopXuPEi8clcmhuzg4d9XjoD7qOf3N7vG52v3rAEM/cCiun1zdqEW4qRQyzz40jKEUNwhvbMivZuaWT+371CirIB6c17F3W6GIy2YpzK0Jl7DBMnffnwVfIrYjR1FAhdNSPEVIM58X9iBhyyPTD6EswFZHaY2VgHEEpmlfVJLyx67WbOTRygV8cHPG4YSXwss58MeIWJujXQrGD0+uLywjJGtL6Ek1t3h979WU6dyUutuGVcFMpMllAxSq3wjiCpahhXvztV6+nvaWB+4LKND53AibO+GjcMZsXLyyc+1Qe2cmtOG/5c3yvsXyaJoOKlf1CxyvhplI4CXtxcYx45AhE5BYR2Ssi+0XkE0U+/5iIDNni9Y+LyGWuz2ZF5Bn756HF3w2VGoQ3Whrr+Z1Xb+KHL5xgdDyA5BLnJvRj6A+28Mb6+Bi3ldNz+S0d/hw/TiWYnZh5v2wD4hU55PdoMYa5FTU7AhGpB74EvB3oB3aJyOIr/DQwYIvXPwD8X67PLiqlrrN/biNKuIU3quCOnZuZmp3ju0EsGvsVMeQmTotgfk6TQbymykb2+xMz7yaT1Zncccit8FK4qRj1DdAdL2U/L0YEO4H9SqmDSqkp4D7gdvcOSqknlFLOysmTaJH66FNjRm3funau713NfbsP+79onM9BW8b7mHk3mX7Ivxh94Y3ZGd1OPx3Byh5o7Y7HAqmfEUMOmRhNh+SHvRNuKkWcRkh44wg2Au6J8CP2tlJ8GPhH1/sWERkUkSdF5J2lviQid9n7DebzAdW88UB4Y9cNvey3zjP4ss8KX35GQTg4uRVRF944/RLMTvrbA4b43OxOzHy3xzHzbgojpBiMGIO6V84eiU1uRaCLxSLyb4EB4POuzZcppQaA9wNfFJErin1XKXW3UmpAKTXQ0+OBulI5NLbo4WMNxv2Oa9ezsrmBe/3MNJ6b83cx0CEuvb4gpslgvuZQ1COH8sPQtU0nO/nFqvV6PSbqtuG1cFMpCiG18Qgx9sIRHAU2u95vsrctQERuBj4F3KaUKqTcKqWO2r8PAj8FrvegTd5Ro/BGa1MDt1+3gUeeP87YxcoXncti7BWYHvdfQLzbya2I+M1u5QCZb69fZPp0rsJYyBoUy+F1Ke5ixCX73LmXA+s0xWCEhDeOYDewXUS2ikgTcAewIPpHRK4H/j+0E7Bc29eISLP9uht4PRCtK+eB8Maunb1MTM/x4DM+JWP5HQXh0LwyHsIbhZj5Vn/PE4cF46kLOrnJb9uAeORW+B1d59DRC41t0bYNFzU7AqXUDPAR4FEgB9yvlNojIp8REScK6PPASuDbi8JEs8CgiDwLPAF8TikVradMJluz8MbVGzu4emM7/+OXPmUa+5k1uphMf/SHu36U4i5GIYQ0wjf7KZ9j5t1k+nUuy/mT/p+rWqxhnfy2+rLl962Fujo9Cot6p8nGk2VzpdQjwCOLtv2p6/XNJb73r8CrvGiDb7h7feuvrfowu3b28qnvvcCzR8a4bvNqb9rmYOWgYzO0tHt73GJksrD/MR2b3tDk//kqZWZSh0tmb/X/XCtW61yFKDuCoEaLsHA6ZNU6/89XDX4IN5Uik4UXf+T/eTzAZBYvR+flUNdYs2e/7doNrGis96c8dRBREA6O8MZodbkVvjOyH9RsgNcj4rkVVg7qm2HNVv/PFYdggqBGi6DXTMYtGA9Bm6RCjCNYDo+EN1a1NHLrtet56NljnJ+c8ahx6Jj5U3uDmRaC6JeaCCpiyKGnz9atiGhuhZXzP2beoa1b6z9E1TbGT/kj3FSKODhGG+MIysGjePFdO3u5MDXLQ88c86BRNqMHYXYquF5OV8SFN6whqGvwvs58KTL9Omdh9KVgzlcpQY4WIdq5FYVOQlCdphgEE9gYR1AOmawnwhvXbV5N37pV3Lfbw+khv1TJStHYEm3hDStnx8wHtH4R5RLME2M6qSlQR9BflbJfIAQVOuqwah20rI7uCMmFcQTl4JHwhohwxw2bee7IGC8cHfOgYczHzDv6CUEQ6V6fzzWGFtNzJSDRvB5WwA8+0Nc+qrkVfgo3FUMkHlF2GEdQHh7Oi7/r+k00N9Tx778+yOf+cZg9x8ZqCym1hqBzKzSuqLltZdOTrTm3whemxuH0y/NlgIOgqQ3WRDS3wu9S3MWIcglmZ6HYD+GmUjghpFHOrcA4gvJYvQUaVnji2TtaG/m7DwywY+0q/u7nB/ntv/1nbvp//om/fuxF9p08V/kBg4yCcHCEN6ImTeh3nflSOKUmooaV00lNHb3BnTOqwQR+CzeVItOvp+jOHQ/2vBUSQChBAvA4OeSNO3p4444eRsen+OELJ3j42WP815/s428f30ffulXceu0G3nHNei7rWkZNamZSl8juv33p/bzGvQi24bpgz70UQcbMu8lkYd+Popdb4ZSWCCJm3qGlA9o3Rc8xnjuuH8ihdJrQ/4v2DcGeuwKMIyiXTD/s/7Gnh+xsa+L9r+3l/a/txTo7wSPPH+f7zx3n84/u5fOP7uWaTR2845r1/PY1G9i4usjUz6l9wcbMFxpuC29EbYHUGtIx850BxMy7cXIrRvbD2oAfNEth5WDH24I/byYbQdsIMPvejXuqbFvRvNpIYBxBuWSy8Mw3YXwE2rq8P3x7Cx98/VY++PqtHD1zkR88d4zvP3ecv3xkmL98ZJjXXLaGW69Zz29ds57Mqhb9pbB6wFEV3rByevG2rj7Y87p7fVFxBOOnYDwfvG2Avh4v/UznuASRv1AOQeeXOLR1wcq1Nech+U1E/ksxoMcVJtj2Bl9PtXH1Cu564xXc9cYrOHRqnB88f5yHnz3Gnz08xJ9/f4jf2NrFrddu4F2jz7OirgE6i1bu9pdMH7zyZPDnXYr8MGzx939TlK5t0cutCKsHDPphOzupdSH81ECoBCsA4aZSxEDP2SwWl0tIWYJbutv4/Tdv44d/9EYe++M38gc3bufk2Qn+8/ee519+8XOO1G/igWctzk74VOK6FJmsDhGMivDGxTNw9mjwPT7Qdf67tkXTEYQ1IoBoPfzCWCh2cEJIo5hbYWMcQbm0b4DmcIU3tq9dxcfesoPH/+Ob+MFH38DAihMMz23if//2swz8xY/5918f5MFnjjLuZQmLUhRyKyIy5A06WWgxUas5ZA3pZKYwir91Ryy3Ym5O20eYtjF9Ac68HM75y8BMDZWLSGQSqUSEq7obYPIYN735Q3xv62/y/eeO8/3njvHY0ElaGuu4qW8tt167nn9zZYaWRh/mzN0jpM07vT9+pQShy7sUmX4YelDX//dbB6EcwoiZd2hq1Qv2EbhXAFu46UKw+RRu3FF2QQcylIlxBJWQ6YM9/1PHJIdxg7mxe8CS6ef63jVc37uGT/1WlsGXT/Pws8d45Pnj/OD547Q11fPWq9Zx67XrecO2HpoaPBoERk14w8pB00pdjjsMnNyKU3thQ8gie0rp6/Gq94TXhijlVoQ5TQbzWf/5HPT9VjhtWAbjCCoh0w9PfQ3OnYD2gNLUS1EkCqKuTti5tZOdWzv59K39PHlwlIefPcYP95zge08fpWNFI7dctY53XLue113eRUN9DU6hrk4beFSmQ6ycXpQLy0G7R0hhO4Jzx2FyLLzREej/xd5/1Lkufmoll0NQqmSlaGnXHZSoOMYieOIIROQW4G+AeuDLSqnPLfq8Gfg68BpgBHifUuqQ/dkngQ8Ds8BHlVKPetEmX3AXGIuCI2ho0ZKMRWior+MN27t5w/Zu/uKdV/PP+/N8/1k9SvjW4GG62pq45ep1bO1uo7Wpgbbmev27qZ7WZv17RVM9bU0NtDbX01Rfhyx+yGb6dSJVFLBycOXbwzv/mq06hyEKjjHsaTLn3GpW57qsuzq8dkCwwk2liMi0cilqdgQiUg98CXgLcATYLSIPLZKc/DBwWim1TUTuAP4KeJ+I9KM1jq8CNgA/FpEdSqloFnd3z/VdcWO4bbGGyo6Zb2qo48a+tdzYt5aJ6Vl+ujfPw88d4zu/PsLEdHmRDA11QmtTPW3NDQUH8d6pFXxg3OLj//AErOjS2y9xKNqRtDbq7zrHaG3S+9XXedCDP5+HC6fCG/qDjpfv2RGNm70QOhqmI3DdK1FwBGE6RdDnP/hTmJ3WGicRw4sRwU5gv1LqIICI3AfczkIR+tuBP7NfPwD8v6K7l7cD9ymlJoGXRGS/fbxfeNAu74mS8IY1DJe/qeKvtTTWc8vV67jl6nXMzinGp2a4MDnLhakZLkzNMj5p/7a3j9vbL0zNMG7vNz41y4XJGV45twWAyWN7eHI2W9h3dq78AlvNDXW0NtXbDkGoE6gTQezfoGeh6kQQ7G32Z3UCgnDdzLP8FfCpf53lucF/RkQvqOv9548nrmPU1envOtv1sbh01IPevuB9Ud8lfPjcWnbkn+Pj9wyW3K/YV4vvV72D/ID1c66q7+Tj3z0EHKr6OLXMstWraf6Gen70xBM89Nzl1R+oRurULF88uZcnJq/ie994KrR27Dy7gg/NTvHn9zzMiaba9JI/fetVrOto8ahlGi8cwUbAXXP2CPDaUvsopWZEZAzosrc/uei7G4udRETuAu4C6O0NsIjWYnr6wu/1XTwN547V3MuprxPaWxppb6myh3J2M/z1p/jim5th500AKKWYnJnj4tS8Eyk4l0nHqSx0LBemZplTijmlv68UzCmFwv5tb59TFLZh7zOnFNefOQ7jcHH1DnrqmwvfmbMrPs4pxdwcKPQxZucU07PKdQ7APv5i9B6u98X2sbfl5jbxutkfc+b0Kcbl0jpR5VaZrbVQZdeFAxyQXl46NV71MRb/3dVwuG4jHef2cWDufM3HqpbNs0doZJpnJtdzIB9eO2R2LR8C6k4Nc6CxtsoEUzPe5yPEZrFYKXU3cDfAwMBAeDVdM/3w9Dd0bHKQxbzcOOnqYQ79wRbe6FgwQhIRWhrraWmsZ01bQAXYHv4HuLCGv/5f3xpuNNfecbj3azzwO2ugd3FfKCDm5uD/PAavvpNH3/7GcNrg8O3XsOXY0/zoDysfuXrGnv8J34Y/+Xfv4k/CLJA4vRM++zH+jxsE3hzi9SiBF0+yo4A7Zm+Tva3oPiLSAHSgF43L+W60yGRhejxc4Y0oLAbCvPBG2COkMGPm3UQho/bMy3bMfMi2Afp/cvqQ1okIizCEm4rRuEIXa4zCtHIRvHAEu4HtIrJVRJrQi78PLdrnIeBO+/V7gJ8oPU5+CLhDRJpFZCuwHfiVB23yjyjokFo5aFoFHZvCa4ODk1EblvCGEzMfhQdfx2adyxC2bUC4C+cOhSi7ELPPrSH9AA5SuKkUEY4cqtkRKKVmgI8AjwI54H6l1B4R+YyI3Gbv9hWgy14M/hjwCfu7e4D70QvLPwR+P7IRQw5REN7ID2ujCrsHDC7hjRPhnP/sUZg8Gw1HUFen15DCLMFciJkPuQcMrk5TiI7AuVeiQKYfRg/A9ETYLbkET9YIlFKPAI8s2vanrtcTwHtLfPezwGe9aEcgtHRA+8bwPLtScHIPZG8N5/yLWSC8EUJuRZR6wKCvx4s/DO/8Vk5nfYcZM++wZovOdQmr0zQ9EY5wUykyfaDmYGQfrHtV2K1ZgCk6Vw1hDvHG83BxNDq9nLA1asMst1yMTFb/j87nwzl/fji8mjqLqauH7hBzK0ZCEm4qRRSmlUtgHEE1ZLJw6kUtvBE0UVkodmjr0nXew3QEK9dBa2c451+MO/s8aGantV1GxTYg3GCCqI0WO6+AusZILhgbR1ANmf554Y2giZpxQ7glmMOsM1+MMHt9owdhdip6tnHumNaLCBprSD94wxBuKkZDkxbqMSOChBBmmKCVg9YuneEcFTL9kN8bvPDG3Kw+b5QefCvXwoo14dzsURstQri6FdawFgxqCCifpRyiplthYxxBNYQpvBGVmHk3hdyKV4I97+lDMHMxWg++MHMrrBxInZ6XjwqhdpoiNloEvaZ25hWYDC/LuRjGEVRDU6uOiAjauJ2Y+agsjDqEJOM5r0oWtZvdLkMSdG6FldNVUKMQM+/QsUnnvARtG5PndXJdlEaL4FpD2htuOxZhHEG1ZPqDj48eOwJT56L54IPgHWOUYubdZLJaD+DssWDPG5XEOjciOoop8E6C/aCN2vWIQvZ5EYwjqJZMFkb2a+GNoIjiQjGEJ7xh5WB1LzSvCva8yxHGgvH0hE5WipptQDjz4lFcLwE7t2JF5BaMjSOoFrfwRlA4IYlRiRN3k8kGP0Jy1kuiRhghpKde1MlKUXvwgf4fXRgJNrciP7ykcFNo1NVHS9nPxjiCagljXtzKwaoNOiolamSyWq83qNyK2WnthKP44Gvt1LkNQdsGRNsxBvnwq0C4KXAy2XDrLxXBOIJq6doOdQ3BG3cURwOgoyFmp3QsexCMHIC56fBLcZci0xesbeRzOma+KyIx827CyD6P6mgR7NyK43BhNOyWFDCOoFoamnSMclCePYox826C7vVFdQ7YIejcCiunk5UiKIPIygys6AzONi6e1g/aKNsGRGpUYBxBLQS5CHb6EMxMRNe4e+zciqCMO4ox824yWa0LcOblYM4XxZh5Bye3IjDbcMKKTaepXIwjqIUghTcKc8ARvdmDFt6whnTpgEZvtVs9I8jIoclzOkkpqrYB84Uag8itKIQVR3QatX0jNLdHKnLIOIJacAwtiJ5O1KpsFiPIqqxRjJl34+Q2BOEYCzHzEe0Bg51bcVbrR/hNlISbiiFiJx2aqaFkEGSvzxqC1ZdB06Wi6JEhk9WLuH4Lb0xf1AX/ouwImldpXYBAbCMmnQQI7npERbipFGEr+y2iJkcgIp0i8piI7LN/XxLXKCLXicgvRGSPiDwnIu9zffY1EXlJRJ6xf66rpT2B07kV6psDNO4I9/hgPrdixOfciijHzLsJaoRk5XSSUtRi5t0ElX2uVLTXSxwy/VpX5LwVdkuA2kcEnwAeV0ptBx633y/mAvABpdRVwC3AF0Vktevz/6SUus7+eabG9gRLITnE55t9Zko/XONg3OD/kDfKMfNuMln9f5ud9vc8UY6Zd2jthFXr/beNgnBTDGwDIrNgXKsjuB24x359D/DOxTsopV5USu2zXx8DLCBCNZRrJIhKk6MHYG4m+sYdlPCGNQT1TXpxOspk+oPJrYjDaBGCibIrhBVHeJoMIqdWVqsjWKuUOm6/PgGsXWpnEdkJNAEHXJs/a08ZfUFEmpf47l0iMigig/l8SDKAxcj02cIbp/07R9Rj5h2CEt6wcjpsNIox826C6PVdGIXzJ6JvG+DKrZj17xxxGS2u7NG6ImEo2RVhWUcgIj8WkReK/CxQhFZKKaDkyoeIrAf+AfiQUsrJsvkk0AfcAHQCHy/1faXU3UqpAaXUQE9PhAYUQUyHWDmQev2QjTo9AWTUWsPRXhh16N6hcx38dIxRLcVdjJ4+rR9x+pB/57CGoifcVIowZTwXsawjUErdrJS6usjPg8BJ+wHvPOiLrnyISDvwA+BTSqknXcc+rjSTwN8DO734owIliAJjVk6XDmgoOWCKDpl+nUTll/DGxFktgBOHB19ji/+5FXEZLUIw0yFRFG4qRZC5FctQ69TQQ8Cd9us7gQcX7yAiTcD3gK8rpR5Y9JnjRAS9vvBCje0Jno7N0LQyAOOOwY0O8+085ZPwRhxi5t34XZXVyunkpPaN/p3DK5zcCr86TUrpax2ne2XqPIwdDrslNTuCzwFvEZF9wM32e0RkQES+bO/zu8AbgQ8WCRP9pog8DzwPdAP/pcb2BI+Iv2GC0xf1YmOcHnzg3/WIUw8Y9P9t1MfcijjEzDs0r9S5MH7ZhiPcFIdpQwinGF8JGmr5slJqBLipyPZB4Pfs198AvlHi+zfWcv7I0NMHwz/QPRKvb8j8XkDFx7jXbNF14H1zBDlobNUPlDjQ06dzHk69COuv8fbYTsx89jZvj+snfnaa4rJQ7OBENlk52PG2UJtiMou9wEkOGfchmiluxu238EYhZj4mpuvnvPj5kzpaLS62AbZuxYs6N8Zr4hI66rBijdYXicCIICZ3U8TxczokLjHzbvyMhsgPx+vB1+VjbkXUCxEWI9Ovc2JGDyy/b6VEWbipFGHIeBbBOAIv8LPXlx+G7iuhvqZZvGBxhDe8zq0YH9G94Dg9+OobdRipH4UJ4zZaBH87TfkYBVU4ZLL+51aUgXEEXuCn8EacIoYc/MqtyMewBwz+9fqsIWjt1slJcaFru86J8doRFISbYmgbs5Mw+lKozTCOwAv8ihyaOKtDy+Iy5+ngV4GxQpXNuN3sfVovYPKct8eNYyfBr9yKqAs3lSKIPKQyMI7AKxxBai+TQwpZozEa+oOuA9+0ynvHaA1Bcwe0b/D2uH5TkCb0MLdCKXu9JGYPPvCn0xS3sGKHHlfkUIgYR+AVfghvxNW4/RohOclCcYiZd+NHzaGxwzoZKW62AXZuxUGdI+MVcdBkKEZTmw65DnnB2DgCr/BjwdgahsY2LXASN7wW3ohLnflirN6i9QK8XDOJ40KxQyYLKB1G6hVWTj9QoyzcVIoI1BwyjsAr/JgXt4b0/HJcYubdeJ1bce4ETJyJ54Ovrk7/H722DYhfDxh86jTl4rd25NDTByP7/cmtKJMYPmEiSmsnrFxnjNsh47FjjFuy0GJ6PJ4qK8TMr/bumEHRebnOjfHKNuIi3FQKJ7diZH9oTTCOwEu8DBMcPwXjVryNG7x7+MV5KgT0//H8Ca0f4AVxnSYDnRPTvcM72xjZHw/hplJEQK3MOAIvyfRD/kVvkkPimDXqps0W3vDKuPM5fcy2bm+OFzReOsa5WW1ncbUN8DaYIK5BFQ7dPuVWVIBxBF6SyXonvBHX0FEHEXsRzKMF0jjGzLvxMl589CWdhBRX2wB9PcYO61yZWskPx0e4qRgNzdC1zTiCxOBlr88agpbVsGpd7ccKC6+EN+bm7NDRGD/42jfoHAivbANi7hg9zK2Ik3BTKbwOJqgQ4wi8xBHe8ORmj1Gd+VL09On68GNHajvO2CswPR7PCBkHEftm98g2YN7e4oiXUXZxXi9xyPTrmYSpC6GcviZHICKdIvKYiOyzfxct+ycisy5Rmodc27eKyC9FZL+IfMtWM4svzSthdW/txh3nmHk3Xo2Q4r5Q7OBVboU1FN+YeYfVl2ldiVptY+qCnipLgm2g/FP2W4ZaRwSfAB5XSm0HHrffF+OiUuo6+8etovFXwBeUUtuA08CHa2xP+GT6a680ee44TIwlwLg96vUVHEGMRwRg51ac1hVUayFupbiLUVenRwW12sYpW7jJdJpqolZHcDtwj/36HrTucFnYOsU3Ao6OcUXfjyxeCG/EPWLIwRHeqNUxWjlo3wQtHd60Kyy8KME8M6nDJeNuG+BNp8mKeVCFw5qtUN8c2jpBrY5grVLquP36BLC2xH4tIjIoIk+KyDvtbV3AGaXUjP3+CFBSgVtE7rKPMZjP+6AE5hVeCG/EtcpmMbzIrYh7xJCDF72+uMfMu8lk9ehofKT6YzjCTWu2eteuMCjkVvigW1EGyzoCEfmxiLxQ5Od2935KKQWUmvy8TCk1ALwf+KKIXFFpQ5VSdyulBpRSAz09Ea6/7sUimJWDtgy0dXnTpjCpVXhjdkYP/+M+LQQ6B6K1u3bbgHgvnDs4/9NaQmqtXPyEm0rhp57zMizrCJRSNyulri7y8yBwUkTWA9i/rRLHOGr/Pgj8FLgeGAFWi4jzH9wEeFi6MyS6d4DU1fYPTcJCsUMmq+vEV5tbMXoQZqeS0QOG2m92ayjeMfNuvBghJWW0CPrvOHtErw8GTK1TQw8Bd9qv7wQeXLyDiKwRkWb7dTfwemDIHkE8Abxnqe/HjsYW6LyieuOem0vGYqBDrenzcVUlK4UzL15t5JA1rJOP4hwz77BqvV73qdY2Jsb0gzNJtgGhTA/V6gg+B7xFRPYBN9vvEZEBEfmyvU8WGBSRZ9EP/s8ppZz//MeBj4nIfvSawVdqbE80qKXXN/YKTF9IjnEXpsqqNG4rB4ge/ieBTFbrCIwdru77SRot1pp97iSjmU5TzdQ0saaUGgFuKrJ9EPg9+/W/Aq8q8f2DwM5a2hBJMv2Qe1gLbzSuqOy7SYmZd6hVeMMags6t0NTqabNCwz0dsrpCnYmpcT3Fdu0uz5sVGpksvPBdPUKqNHky7hVpF9OxWeuP1BpJVQUms9gPMn2Aqi59vlBnPiE9YKitBHOcS3EXo5B9XoVjzDsx8wl58IH+306c0XoTlWLl4ivcVAw/dCvKPXXgZ0wDtSyCWTndM2hp97ZNYZLJ6nrxleZWzEzCyIHkTIWA1g9o31i9bUByRotQ23RInIWbShFS5FCCrmCEcIQ3qgmLS1IUhEO1whun9oGaTeD1qDK3Ip/TSUdxj5l3U0uSXVLvlfE8nA82V8o4Aj+ob6xOeGN2RmclJ864qyzBnMQeMNi5FVXoVlg56NmRjJh5h7ZunTNTqW2Mn9IPzCTaBnhTrrwCjCPwi2qGeEmLmXeoVnjDGoK6Bh0umSQy/VpPYPSlyr5n5ZJnG1DdvZKkxDo3znpYwCGkxhH4RU9f5cIbcRYkX4qGZl0vvpqbvWsbNMS7KO0lVJN9fvEMnD2aPNsA2xEM6xyacknqaHHVOq1DEvCCsXEEflEQ3qjAszsx80mKGHKoZl48STHzbnquBKQyxxh3xbqlyGS13sTYK+V/JwnCTcUo5FaYqaFkUM0imDWkF5orzT2IA5l+PRVSrvDG1DiceTmZD75qciuSUpG2GNVE2TnTZHEWbiqFV8p+FWAcgV9UI7yRH07mjQ4u4Y0Xy9u/0ANO6vWosASzlYOmlTq0OGkUpsrKvFeU0oupibWNLEyOwdljgZ3SOAK/qFR4Y3rCjplPYA8YKu/1JXUO2CGT1eG0M5Pl7W8NaXtKUsy8Q0u7dnDl2kZBuCnBjgACjRxKoFVFiEqiIUacmPkELgaCLbzRVL5jtHLQ0KKnUJJIJltZboWVS65tgN1pKreT4JSWSKgj6Kkht6JKjCPwk0wWxq3yhDeS3gOub9CF4yq52bt3QF29v+0Ki0rWkM7n4cKp5NoG2Mp+e3UuzXIkSbipGG1dsHKtcQSJoZIhnpWDukZdwjqpZLLlz4tbCSrFXYyu7TpHopybPWmluIuR6dc5NKfLyK2whvWDMgnCTaXwQtmvAowj8JNK5sWtnE68SlrMvJtMtrzcioun4dyxZD/4Gpp0jkS5tgHJdoyV1BxKalixmx5H2a+C3IoaMI7AT1ath+YyhTfSYNzl5lYkRZB8Ocrt9VlDsGKN7gUnlXJzKxzhpqROCzlkslqX5MzLgZzOOAI/ESlvwXjyvP6HJ964y8yoTVqd+VL0ZLW+wNT40vs5pbiTGDPv0LhC604sZxtnXk6WcFMpvJDxrICaHIGIdIrIYyKyz/69psg+bxaRZ1w/EyLyTvuzr4nIS67PrqulPZHE6fUtlRxSUFpKuHF39Or68csZd5Jj5t04uRVL6VYolcwqm8UoJ6M2DdNkUJtuRRXUOiL4BPC4Umo78Lj9fgFKqSeUUtcppa4DbgQuAD9y7fKfnM+VUs/U2J7okenXMc9LCW+kYTEQXMIby9zsTmJdknvAUN5U2dljMHk2+bYBdm7FgaVzK5x7JYllWNy0tOuOUxxGBMDtwD3263uAdy6z/3uAf1RKlVlnIAGUswhm5aBhRXJj5t0sN1WmFJzck44HX+dWrS+wnG1A8nvAoP/nalbrUJTCyukHZJKEm0pRTqfJI2p1BGuVUsft1yeA5Vaz7gDuXbTtsyLynIh8QUSaS31RRO4SkUERGczngxVtqIly4sWtId3DSWrMvJtMv51bcar45+N5uDiajgdfXb3+vy9nG5AOx1jOvHjSE+vcOMp+s9O+n2pZRyAiPxaRF4r83O7eTymlgJIT4SKyHi1i/6hr8yeBPuAGoBP4eKnvK6XuVkoNKKUGenp6lmt2dGjrhraeMow7BTc6LF9XJqmluEux3AjJyuloodbO4NoUFp1X2LkVJUZIs9PJFG4qhZNbMXrQ91Mt6wiUUjcrpa4u8vMgcNJ+wDsPemuJQ/0u8D2lVMG9KaWOK80k8PfAztr+nIiyVJjgxdO6dkqajBuWcAQpmgoB/X8/e1TrDRQjDWHFDg1NOtGulG0kVbipFLXoOVdIrVNDDwF32q/vBB5cYt9dLJoWcjkRQa8vvFBje6JJpr90ckhaYuYdHOGNUtnWVg5WdMLKTKDNCo3CgnGRyKG5Ob09LbYBdvb5cp2ElDjG7h0gdYGsE9TqCD4HvEVE9gE32+8RkQER+bKzk4hsATYD/7To+98UkeeB54Fu4L/U2J5ospTwRprmgGF54Y0k15kvxlK9vjOHYOZiemwD9P++VG6FldMPxu4dgTcrFBpX6GKNAYwIalLBVkqNADcV2T4I/J7r/SFgY5H9bqzl/LHBXU1wcWSQlYPmdmi/5PIkl0wfvPAdHSHkfuA7MfPXvi+8tgVNx2adM1HMMSa9uFoxnIXg/DBsfM3Cz6wh/WBMonBTKRwZT58xmcVBsFRGrZXTC6Np6QGDK7fi+MLtY0dg6ly6esAipXUrCgvnCY+Zd7PUGlKagiocMv0wekDrlfiIcQRB0NIB7ZsuNW6l0rUY6FBqOiRtC8UOpSKH0hQz77Bmi9ahWHw9pif0AzGNtqHmylf2qxLjCIKi2BAvTTHzbgpTZYuuRyFrNCWhow6Zfq03cH5RfoyVYOnSUpTKrRjZpx+IabseAdUcMo4gKIoJb6RtodihlPCGldMVW9MQM++mmG5F2mLm3RQLJkjraLHrCq1T4vOCsXEEQZHJXpockrZwODfF5sUdXd60USz7fOQAzE2n1zbOHdM5Ng7WkH4gdiVYuKkY9Y1ap6RcQacqMY4gKIrNi1tD0NqlM4/TRqZfG7eTWzE3m76YeYeVa7XewGLbgHQ6gsJ0iOvh5wg31TeG06YwCUCtzDiCoOguIryRtph5N4uFN04fgpmJdD74iuVWpC1m3k2pTlMabQP0333mFZg859spjCMIiqZWXW3SmQdWKp2LgQ6LSzCndQ7YwQkmcHQr8jnovDxdMfMOHZugadW8bUye1w/C1N8rS+hW1IhxBEHi7vWlMWbezWLhjULyVIpi5t1ksjA5pvUHIJ0x8w6Llf0Kwk0p7ST0LJGH5BHGEQRJT59eBJyeSGfWqJuWdp1V61wHawhW90LzynDbFRbu7PPpizqoIK22AToJ8+Se+VwbSGcgAdi5FSt8zTA2jiBIHOGNkX3p0eVdCnevz1kvSSvuefFTL6YzZt5Npl/n2Izn0yXcVIxCboUZESQDd3KIlYNVG3S0SFrJZPVDb/qido5pfvC1dsLKdfO2AcYxgn74pUm4qRTl6DnXgHEEQdK1zRbeyOnFwDQ/+GBeeGPfj2BuJt0PPpgvwWzl0hkz78YdQpofNraRycL5E3Bh1JfDG0cQJI7wxskX7Jj5tDsC++9/4TsL36eVTL9+8J3co8NG0xgz79DWo3NsXv6XdAk3laIcydsaMI4gaDJ98NLP0hsz78YR3njxUf27a3vYLQqXTJ/WHzj083SvHYFdlTWrR4tg7pViZUg8xDiCoMn0aycAxrgd4Y2ZCa1X29gSdovCxZn+MJ0ETSZr7hWH9o1atySKIwIRea+I7BGROREZWGK/W0Rkr4jsF5FPuLZvFZFf2tu/JSJNtbQnFrgNOq3hcG6c65H2Gx0W5lCkfU4c5m0ibcJNxVicW+ExtY4IXgB+B/hZqR1EpB74EvB2oB/YJSKOlf8V8AWl1DbgNPDhGtsTfZwbfM0WaGoLtSmRwLke5sEHzat0LgUYxwgu28imswzLYpxCjU72uYfU5AiUUjml1HJ5zzuB/Uqpg0qpKeA+4HZbsP5G4AF7v3vQAvbJxhHeMA8+jRkRLCTTr2PmV28JuyXh46yTGNvQZPp1RdbzJz0/dE2axWWyETjsen8EeC3QBZxRSs24tpcc/4nIXcBdAL29vf60NAjq6uFtf5neUgqL2f4WeN1HYNsl0tfp5Df/ALa/FerM8h0r1sDNfw5XvDnslkSDTQNw1bvm1008ZFlHICI/BtYV+ehTSqkHPW9RCZRSdwN3AwwMDHg/NgqSG5I/A1Y2zavgbZ8NuxXRYcsb9I9B84Y/CrsF0WHTALz3a74cellHoJS6ucZzHAU2u95vsreNAKtFpMEeFTjbDQaDwRAgQYw/dwPb7QihJuAO4CGllAKeAN5j73cnENgIw2AwGAyaWsNH3yUiR4DXAT8QkUft7RtE5BEAu7f/EeBRIAfcr5TaYx/i48DHRGQ/es3gK7W0x2AwGAyVI8qHUCS/GRgYUIODg2E3w2AwGGKFiDyllLok58uEJhgMBkPKMY7AYDAYUo5xBAaDwZByjCMwGAyGlBPLxWIRyQMvV/n1buCUh82JO+Z6zGOuxULM9VhIEq7HZUqpnsUbY+kIakFEBoutmqcVcz3mMddiIeZ6LCTJ18NMDRkMBkPKMY7AYDAYUk4aHcHdYTcgYpjrMY+5Fgsx12Mhib0eqVsjMBgMBsNC0jgiMBgMBoML4wgMBoMh5aTKEYjILSKyV0T2i8gnwm5PWIjIZhF5QkSGRGSPiPxh2G2KAiJSLyJPi8j3w25L2IjIahF5QESGRSQnIq8Lu01hISJ/bN8nL4jIvSLSEnabvCY1jkBE6oEvAW8H+oFdIpJW4eAZ4D8qpfqB3wB+P8XXws0fokulG+BvgB8qpfqAa0npdRGRjcBHgQGl1NVAPVpTJVGkxhEAO4H9SqmDSqkp4D7g9pDbFApKqeNKqV/br8+hb/KSetFpQEQ2Ab8NfDnstoSNiHQAb8TWB1FKTSmlzoTaqHBpAFaISAPQChwLuT2ekyZHsBE47Hp/hJQ//ABEZAtwPfDLkJsSNl8E/gSYC7kdUWArkAf+3p4q+7KItIXdqDBQSh0F/m/gFeA4MKaU+lG4rfKeNDkCwyJEZCXwHeCPlFJnw25PWIjIOwBLKfVU2G2JCA3Aq4H/ppS6HhgHUrmmJiJr0DMHW4ENQJuI/NtwW+U9aXIER4HNrveb7G2pREQa0U7gm0qp74bdnpB5PXCbiBxCTxneKCLfCLdJoXIEOKKUckaJD6AdQxq5GXhJKZVXSk0D3wV+M+Q2eU6aHMFuYLuIbBWRJvSCz0MhtykURETQ8785pdRfh92esFFKfVIptUkptQVtFz9RSiWu11cuSqkTwGERudLedBMwFGKTwuQV4DdEpNW+b24igQvnDWE3ICiUUjMi8hHgUfTK/1eVUntCblZYvB74d8DzIvKMve0/K6UeCa9JhojxB8A37U7TQeBDIbcnFJRSvxSRB4Bfo6PtniaBpSZMiQmDwWBIOWmaGjIYDAZDEYwjMBgMhpRjHIHBYDCkHOMIDAaDIeUYR2AwGAwpxzgCg8FgSDnGERgMBkPK+f8BwaSPXwFUoZ4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Run the 1D model to see where the steps break down..\n",
    "import numpy as np\n",
    "import scipy.integrate\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.linalg\n",
    "\n",
    "def degradation_model1D( t, y, params ):\n",
    "    dydt = -params[0] * y;\n",
    "    return dydt\n",
    "\n",
    "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",
    "    y = np.zeros( (1,t.size) )\n",
    "    y[0,0] = y0;\n",
    "    for i in range(1,t.size):\n",
    "        h = t[i]-t[i-1];\n",
    "        yprime = fn(t[i-1],y[0,i-1],params);\n",
    "        y[0,i] = y[0,i-1] + h*yprime;\n",
    "    \n",
    "    return (t,y)\n",
    "\n",
    "dt = 1\n",
    "y0 = 1\n",
    "params = (2,)\n",
    "t = np.arange(0,10,dt)\n",
    "tru_y = np.exp( -params[0]*t)\n",
    "(t,y) = forwardEuler1D( degradation_model1D, t, y0=y0, params=params)\n",
    "\n",
    "plt.plot(t, tru_y)\n",
    "plt.plot(t, y.T);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rewrite the forward Euler method to work for multi-state models\n",
    "\n",
    "Rewrite the method so that y0 can have more than one value and fn can return a vector of more than one value.\n",
    "\n",
    "The return value should be NY x NT (number of states by number of time steps) big, so it gets filled in one column at a time. Adjust the 1D code to account for multiple states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forwardEuler( fn, t, y0, params ):\n",
    "    \"\"\" [t,y] = forwardEuler( 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 (NY,NT)\n",
    "    \"\"\"\n",
    "    # Write code here\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, make a 2-state degradation model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def degradation_model( t, y, params ):\n",
    "    # Write code here\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the 2-state degradation model with forward Euler.\n",
    "\n",
    "Use 2 and 2.5 as the degradation rates. Simulate from t=0 to t=4 with a time step size of 0.1.\n",
    "\n",
    "Then plot the true solution and the approximate solution (indicate simulated with dashed lines). And in another plot, show the difference (between the true and approximate solutions) over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determine how the error depends on the step size\n",
    "\n",
    "When we compute the difference between the true and approximate solutions, we have one value at each state for each time step. I would like to summarize this as one number, so that we can plot \"error\" versus step size.\n",
    "\n",
    "Let's explore different measures for the error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine how the error changes with the time step. Plot two figures - one with the simulations and one with the errors. Compute the error using two methods:\n",
    "- Compute the Euclidean distance between the true solution and the approximate solution (2-norm of the difference) at each time point, then take the mean across time.\n",
    "- Compute the infinity norm (max value) of the difference between the true solution and the approximate solution  at each time point, then take the mean across time.\n",
    "\n",
    "Make an array of time steps. For each time step, plot each solution in a separate subfigure, and keep track of the error in an array. Then plot the errors in a second figure, with each error computation in a separate subfigure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# range from smallest timestep that doesn't take way too long to run to largest time step with qualitatively correct answer.\n",
    "dts = np.linspace(1e-4,1e-1,10)\n",
    "err1 = np.zeros( (dts.size,) )\n",
    "err2 = np.zeros( (dts.size,) )\n",
    "params = (2,2.5)\n",
    "y0 = np.array((1,1))\n",
    "\n",
    "for i in range(dts.size):\n",
    "    t = np.arange(0,4,dts[i])\n",
    "    sol = scipy.integrate.solve_ivp(fun=degradation_model, t_span=(0,150), y0=y0, method='LSODA', rtol=1e-8, atol=1e-8, t_eval=t, args=[params] )\n",
    "    (t,y) = forwardEuler( degradation_model, t, y0, params)\n",
    "    plt.subplot( dts.size, 1, i+1 )\n",
    "    plt.plot( t, sol.y.T )\n",
    "    plt.plot( t, y.T, '--')\n",
    "    err1[i] = np.mean(scipy.linalg.norm( sol.y-y, 2, axis=0) )\n",
    "    err2[i] = np.mean(scipy.linalg.norm( sol.y-y, np.inf, axis=0) )\n",
    "    \n",
    "plt.figure()\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot( dts, err1, '.-' );\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot( dts, err2, '.-' );\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus Content\n",
    "\n",
    "Here, we use Forward Euler to solve the Goodwin oscillator (coded up in Lecture 8). This code is hear as an example to follow when working on Project 3. Aside from being more interesting than the degradation model, the most important difference is that we don't have an analytical solution. So, we use solve_ivp's LSODA method with really tight error controls (i.e. we allow only tiny little errors) to generate the reference \"true\" solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here is the model, copied from a notebook from earlier in the course\n",
    "def goodwin_oscillator( t, statevars, params ):\n",
    "    X = statevars[0];\n",
    "    Y = statevars[1];\n",
    "    Z = statevars[2];\n",
    "    alpha_1 = params[0]\n",
    "    K = params[1]\n",
    "    n = params[2]\n",
    "    gamma_1 = params[3]\n",
    "    alpha_2 = params[4]\n",
    "    gamma_2 = params[5]\n",
    "    alpha_3 = params[6]\n",
    "    gamma_3 = params[7]\n",
    "        \n",
    "    dydt = np.zeros( (3,) )\n",
    "    dydt[0] = alpha_1*K**n/(K**n+Z**n) - gamma_1*X\n",
    "    dydt[1] = alpha_2*X - gamma_2*Y;\n",
    "    dydt[2] = alpha_3*Y - gamma_3*Z;\n",
    "\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulate the model for a given time step size.\n",
    "# SRT reran this many times, adjusting the time step.\n",
    "import scipy.integrate\n",
    "\n",
    "alpha_1 = 1\n",
    "K = 0.61\n",
    "n = 11\n",
    "gamma_1 = 0.15\n",
    "alpha_2 = 1\n",
    "gamma_2 = 0.15\n",
    "alpha_3 = 1\n",
    "gamma_3 = 0.25\n",
    "\n",
    "params = (alpha_1,K,n,gamma_1,alpha_2,gamma_2,alpha_3,gamma_3)\n",
    "\n",
    "# Find a point on the limit cycle and start all subsequent simulations from that point.\n",
    "sol = scipy.integrate.solve_ivp(fun=goodwin_oscillator, t_span=(0,300), y0=(.1,.1,.1), method='LSODA', rtol=1e-8, atol=1e-8, args=[params] )\n",
    "y0 = sol.y[:,-1]\n",
    "\n",
    "# Run this with a time step of 0.1 (which is the biggest that looks like it is still qualitatively correct)\n",
    "# Run this with a smaller time step - the smallest you are willing to wait for, as well. \n",
    "t = np.arange(0,150,.1)\n",
    "sol = scipy.integrate.solve_ivp(fun=goodwin_oscillator, t_span=(0,150), y0=y0, method='LSODA', rtol=1e-8, atol=1e-8, t_eval=t, args=[params] )\n",
    "(t,y) = forwardEuler( goodwin_oscillator, t, y0, params)\n",
    "\n",
    "plt.plot( sol.t, sol.y.T );\n",
    "plt.plot( t, y.T, ':' );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine how the error changes with the time step. Since we don't know exactly what the region of absolute stability is, we want to make sure we plot every solution (and not just the error we compute from it). So plot each solution in one figure, but keep track of the error in an array. Then plot the error in a second figure. Here, I compute the error by first computing the Euclidean distance between the true solution and the approximate solution at each time point, then take the mean across time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determine how the error changes with the time step.\n",
    "# range from smallest timestep that doesn't take way too long to run to largest time step with qualitatively correct answer.\n",
    "dts = np.linspace(1e-4,1e-1,10)\n",
    "err = np.zeros( (dts.size,) )\n",
    "\n",
    "for i in range(dts.size):\n",
    "    t = np.arange(0,150,dts[i])\n",
    "    sol = scipy.integrate.solve_ivp(fun=goodwin_oscillator, t_span=(0,150), y0=y0, method='LSODA', rtol=1e-8, atol=1e-8, t_eval=t, args=[params] )\n",
    "    (t,y) = forwardEuler( goodwin_oscillator, t, y0, params)\n",
    "    plt.subplot( dts.size, 1, i+1 )\n",
    "    plt.plot( t, sol.y.T )\n",
    "    plt.plot( t, y.T, '--')\n",
    "    err[i] = np.mean(scipy.linalg.norm( sol.y-y, 2, axis=0) )\n",
    "    \n",
    "plt.figure()\n",
    "plt.plot( dts, err, '.-' );"
   ]
  }
 ],
 "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
}