01_Forward_Euler.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "![Logo TUBAF](https://tu-freiberg.de/sites/default/files/media/freiberger-alumni-netzwerk-6127/wbm_orig_rgb_0.jpg)\n",
    "\n",
    "Exercise material of the MSc-level course **Numerical Methods in Geotechnical Engineering**.\n",
    "Held at Technische Universität Bergakademie Freiberg.\n",
    "\n",
    "Comments to:\n",
    "\n",
    "*Prof. Dr. Thomas Nagel  \n",
    "Chair of Soil Mechanics and Foundation Engineering  \n",
    "Geotechnical Institute  \n",
    "Technische Universität Bergakademie Freiberg.*\n",
    "\n",
    "https://tu-freiberg.de/en/soilmechanics\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 1 - Forward Euler\n",
    "\n",
    "The ODE $\\dot{y} = k (y_\\text{eq} - y)$ is to be integrated numerically.\n",
    "\n",
    "First, we import some libraries for numerical basics and for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#Some plot settings\n",
    "plt.style.use('seaborn-deep')\n",
    "plt.rcParams['lines.linewidth']= 2.0\n",
    "plt.rcParams['lines.color']= 'black'\n",
    "plt.rcParams['legend.frameon']=True\n",
    "plt.rcParams['font.family'] = 'serif'\n",
    "plt.rcParams['legend.fontsize']=14\n",
    "plt.rcParams['font.size'] = 14\n",
    "plt.rcParams['axes.spines.right'] = False\n",
    "plt.rcParams['axes.spines.top'] = False\n",
    "plt.rcParams['axes.spines.left'] = True\n",
    "plt.rcParams['axes.spines.bottom'] = True\n",
    "plt.rcParams['axes.axisbelow'] = True\n",
    "plt.rcParams['figure.figsize'] = (8, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first code the analytical solution of the ODE in order to have something to compare against.\n",
    "\n",
    "$\\displaystyle \\frac{\\text{d}y}{y_\\text{eq}-y} = k\\,\\text{d}t \\quad \\rightarrow \\quad \\int \\limits_0^t k\\,\\text{d}\\bar{t} = \\int \\limits_{y_0}^y \\displaystyle \\frac{\\text{d}\\bar{y}}{y_\\text{eq}-\\bar{y}}$\n",
    "\n",
    "$kt = \\displaystyle \\left. - \\ln (y_\\text{eq}-\\bar{y}) \\right|^y_{y_0}$\n",
    "\n",
    "$\\exp (-kt) = \\displaystyle \\frac{y_\\text{eq}-y}{y_\\text{eq}-y_0}$\n",
    "\n",
    "$y = y_\\text{eq} - (y_\\text{eq} - y_0) \\exp (-kt)$\n",
    "\n",
    "Let's define the equilibrium value and the rate constant as global variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_eq = 1.\n",
    "k = 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def y_analytical(t,y_0):\n",
    "    return y_eq - (y_eq - y_0)*np.exp(-k*t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to the ODE definition $\\dot{y}=f[y(t),t]$ we define the rate function itself:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dy_dt(y):\n",
    "    return k*(y_eq - y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the function and its rate for some exemplary values ($y_0 = 0$) in a time interval $t \\in [0,4]$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "time = np.linspace(0,4,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$y$, $\\\\dot{y}$ / s$^{-1}$')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAF9CAYAAACZG2t7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+A0lEQVR4nO3deZxcVZn/8c9TS+/pdNLZOoQQQiCQECGkEcISI4IgICo6/hzU0ZlBGH84Ekd0xo0ZFfjJzOiAjIPDjBgXGHdQHIZNJ0QRMCskhARIIPvSSWftvaue3x9V3el0eu+qure6vu/Xq15V99xz730ul049de4555q7IyIiIoUrEnQAIiIiEiwlAyIiIgVOyYCIiEiBUzIgIiJS4JQMiIiIFDglAyIiIgUuFnQAQbniiiv8scceCzoMERGRXLKeCgu2ZWDv3r1BhyAiIhIKBZsMiIiISIqSARERkQKnZEBERKTAKRkQEREpcEoGREREClzBDi0UEZGRr62tjW3bttHc3Bx0KFkXjUapqqpi3LhxRCKD+62vZEBEREasbdu2MWrUKKZNm4ZZj0PsRwR3p62tjd27d7Nt2zamTp06qO11m0BEREas5uZmqqurR3QiAGBmFBUVccIJJ9DQ0DDo7UOZDJjZJ8zMzWxhP/UqzWyxmS0zs5VmdqeZqbVDREQ6jfREoKvB3h7o3C7DcQybmU0GPjPA6ouBqLufC5wPXAx8JUuhiYiIjEihSwaAe4A7+qtkZmcC7wH+EcDdW4G7gEVmVpHNAEVEREaSUCUDZvZOoA14fADVLwWagbVdypYBpcBFmY+ud6/ue50HXniI57auzOVhRUREMiI0yYCZlQO3A58a4CbTgd3u7l3KdnVZ19MxbjCz5Wa2vK6ubujBdrPlwHZ+uf4JVu5c239lERGRkAlNMgB8Ffi2u+8cYP1yoKVbWcdyWU8buPt97l7r7rXjx48fYpjHG11SCcDB5kMZ26eIiEiuhCIZMLNzgPOAbw9iswaguFtZx3JjJuIaqKp0MnCgScmAiIgM35QpU/jGN75xTNmaNWsoKSlh3bp1GT9eWIbhXUXqXv9v00NAStLld5nZAeBGd9/QbZtNwAQzsy63Cial3zdmOd5jVJWmk4EWJQMiImH2zk//MpDjPvL1dw2q/vz581m2bNkxZYsWLeL6669n1qxZmQwNCEnLgLt/1d3PcfeF7r4Q+EB61aJ0WfdEAOBJUgnE7C5ltUAT8ExWA+5mdPEoAA42HybpyVweWkRERqDuycDDDz/MqlWr+PKXv5yV44WlZaBfZnYb8G6g1t2b3f0lM3uI1JwEHzGzOHAzcJe7H8llbPFonPKiMhpaGznS0kBlyahcHl5ERAZosL/Qg3L++efz6U9/mvr6esrLy7nlllu49dZbqa6uzsrxQpcMmNldpCYQgtRtglfd/U9I3TooA7pOJfVR4B4zWwZEgaeAW3MX7VFVJZU0tDZyoPmQkgERERmWefPmUVRUxPLly1m1ahWxWIybbropa8cLXTLg7ot6Kb8FuKVb2SHgIzkIq19VJZVsP7SLA82HmMoJQYcjIiJ5rLi4mLlz5/LII4/wve99jwcffJB4PJ6144Wiz8BI0DG88ICGF4qISAbMnz+fb33rW8yfP5+rr746q8dSMpAhVUoGREQkg84++2wikchxQwyzQclAhlRp4iEREcmgH/7wh9x4443Mnj27/8rDFLo+A/lKLQMiIjJcyWSSuro6Fi9ezNq1a/nJT36Sk+MqGcgQJQMiIjJcS5cu5ZJLLmHmzJn8/Oc/Z8yYMTk5rpKBDFEyICIiw7Vw4UKSydxPXqc+AxmiPgMiIpKvlAxkSGXxKMyMwy0NtCcTQYcjIiIyYEoGMiQSiVBZVIHjHGo5HHQ4IiIiA6ZkIIOO3ipQMiAiIvlDyUAGdT7KuPlgwJGIiIgMnJKBDOqckrhJnQhFRCR/KBnIIA0vFBGRfKRkIIM0vFBERPKRkoEMUsuAiIjkIyUDGaRkQEREMmnLli1UVFSwatWqrB5H0xFn0GglAyIikkGTJ09m9erVnHjiiVk9jpKBDFKfARERyaRYLMaMGTOyfhzdJsig8qIyopEoDW1NtCbagg5HRERkQJQMZFDEIowuHgWodUBERPKHkoEM05TEIiIyXFOmTOEb3/jGMWVr1qyhpKSEdevWZfx4SgYy7OiIAk1JLCIiQzN//nyWLVt2TNmiRYu4/vrrmTVrVsaPpw6EGabhhSIi4fX+H388kOP+5P/cO6j68+fP59/+7d86lx9++GFWrVrFT37yk0yHBqhlIOM0vFBERIbr/PPPZ+PGjdTX19PS0sItt9zCrbfeSnV1dVaOp5aBDFPLgIhIeA32F3pQ5s2bR1FREcuXL2fVqlXEYjFuuummrB0vNMmAmS0AFgFjgShQBfynu9/dxzZLeih+1t0/l4UQB+ToY4yVDIiIyNAUFxczd+5cHnnkEb73ve/x4IMPEo/Hs3a80CQDwHXAanf/CoCZnQWsNLON7v7r3jZy94U5im9AOkcT6DHGIiIyDPPnz+fuu+/msssu4+qrr87qscKUDHwT2Nqx4O4vmNkBIPtTL2VQZ5+BFg0tFBGRoTv77LOJRCLHDTHMhtAkA+7eOXDSzCLAXwItwE8DC2oI1GdAREQy4Yc//CE33ngjs2fPzvqxQpMMdDCzLwKfAOqBK919ez/17wbmAgY8B9zm7oEN8i+NlVAUjdPS3kJzWzMl8ZKgQhERkTyTTCapq6tj8eLFrF27NmtDCbsL3dBCd78NqAFuB542swv6qP4C8D/uvgC4gtQthSVm1mMvCzO7wcyWm9nyurq6TIfecQy1DoiIyJAsXbqUmpoaFi9ezM9//nPGjBmTk+OGrmUAwN0deMDMPgB8DVjQS72bu3xuMLNbgNeAq4CHe6h/H3AfQG1trWc+8pTRJZXsadjHgebDTBo1IVuHERGREWbhwoUkk8mcHzc0LQNmVtRD8TpgMDdLXgecgDsdakpiERHJJ6FJBoAVZmbdyiYDPfYZMLMJZvaFbsU1pPoObMlCfAOm2wQiIpJPwpQMjCLVcRAAM5sHvA/4Tnr5RjPbbGY16SplwN+Y2bT0+ghwK/AG8Gjuwj6epiQWEZF8EqY+A58Hrjez64AEUAp8GuiYO7IoXRZNL+8Cvg78l5k1A+XAZuBt7n4kl4F3p8cYi4hIPglNMuDuDwIP9rH+HuCeLsvNwB3pV6h0JAP71WdARCRw7s7xd6FHplT/+8EL022CEWNsaRUA+5sOBBqHiEihi0ajtLW1BR1GzjQ1NQ3pGQZKBrJgbFkVAPWNBwKNQ0Sk0FVVVbF79+5AhuvlkrvT2NjI9u3bmTBh8EPaQ3ObYCSpKq4kYhEOthymLdFGPJq9J02JiEjvxo0bx7Zt29iwYUPQoWRdPB5n4sSJVFZWDnpbJQNZEIlEGFM6mn2N+9nfdJAJFeOCDklEpCBFIhGmTp0adBihp9sEWVJdmppCcl/T/oAjERER6ZuSgSzp6DewT/0GREQk5JQMZElHy0C9WgZERCTklAxkSbVGFIiISJ5QMpAlHXMN7NNcAyIiEnJKBrJkbMdtgkbdJhARkXBTMpAlHbcJ1DIgIiJhp2QgS8aUjMYw9jcfJJFMBB2OiIhIr5QMZEksGmN0ySjcXY8yFhGRUFMykEWdEw+p34CIiISYkoEsGtMxvFD9BkREJMSUDGRRdXp4oZIBEREJMyUDWVRdptsEIiISfkoGskgTD4mISD5QMpBFHS0DmnhIRETCTMlAFlWrZUBERPKAkoEsGtulA2HSk8EGIyIi0gslA1lUFCuioqicRDLBoZYjQYcjIiLSIyUDWdY5vFD9BkREJKSUDGTZ2I5OhOo3ICIiIaVkIMs6OxE2Hgg0DhERkd7Egg6gg5ktABYBY4EoUAX8p7vf3cc2ceBO4C2AA6uARe7ekO14B6qjZWBfk24TiIhIOIWpZeA6YLW7L3T3i4EPAd8ws6v72OZOYC5wHnAuqQTiP7Id6GAc7TNwINA4REREehOalgHgm8DWjgV3f8HMDgAzeqpsZmOATwDvdff2dNk/Ac+b2a3u/lr2Q+5ftVoGREQk5ELTMuDu69z9MICZRczsY0AL8NNeNnkLEAeWdylbBSSAS7MZ62CMVcuAiIiEXJhaBgAwsy+S+sVfD1zp7tt7qTqdVD+BXR0F7t5mZvvS60JhbPoxxvua9uPumFmwAYmIiHQTmpaBDu5+G1AD3A48bWYX9FK1HGhzd+9W3gKU9bSBmd1gZsvNbHldXV3GYu5LWbyU0lgJrYk2Globc3JMERGRwQhdMgDgKQ8AS4Gv9VKtAYjb8T+1i4Eev3Xd/T53r3X32vHjx2cu4H50tA5orgEREQmj0CQDZlbUQ/E6YHYvm2wCDJjYZR8xoBrYmPEAh6G6VJ0IRUQkvEKTDAAreviVPxnorc/A00ArUNulbC6pOQqeynx4Q9fZb0CdCEVEJITClAyMItVxEAAzmwe8D/hOevlGM9tsZjUA7r4f+BbwKTOLpROJW4D/cvdQtgzUq2VARERCKEyjCT4PXG9m15EaHlgKfBq4N72+KF0W7bLN35GaeOh5UiMLVgM35yjeARurKYlFRCTEQpMMuPuDwIN9rL8HuKdbWSvwqSyHNmzVnR0I1TIgIiLhE6bbBCNWxyyEexuUDIiISPgoGciB8WXVANQ17uP4aRFERESCpWQgB8qKSikvKqM10cbB5kNBhyMiInIMJQM5MqE81Tqwp2FfwJGIiIgcS8lAjkwoHwcoGRARkfBRMpAjR1sG9gYciYiIyLGUDOTIeN0mEBGRkFIykCMdtwnq1DIgIiIho2QgRyZUpFsGjqhlQEREwkXJQI5MSM81sLexnmQyGXA0IiIiRykZyJGiWBFVJZUkPEl904GgwxEREemkZCCHxmtEgYiIhJCSgRzSxEMiIhJGSgZySBMPiYhIGCkZyCFNPCQiImGkZCCHOvoM1KllQEREQkTJQA5NqEjfJtBcAyIiEiJKBnJoXNlYzIz6pgO0J9qDDkdERARQMpBTsUiU6tIxOM7exvqgwxEREQGUDOScHlgkIiJho2QgxzSiQEREwkbJQI5p4iEREQkbJQM5pomHREQkbJQM5FjnXANHdJtARETCQclAjk2o0G0CEREJl1jQAXQws6uB/wsUA+XAfuBv3f3FfrZbD+zqVvyQu9+dlUCHaWxJFdFIlIMth2lub6EkVhx0SCIiUuBCkwwAi4FPuvuDAGb2NeA3Znamu+/uY7td7r4wB/FlRCQSYXzZWHYdqWNvQz1TRtcEHZKIiBS4MN0mWNqRCKR9HRgHvD2geLLmaCdC9RsQEZHghSYZcPdruxU1pd9HXDu6Jh4SEZEwCdNtgu7mA83Ar/qpV25m9wMzgCTwFPDP7t6c5fiGrHOuAY0oEBGREAhNy0BXZmbAl4AvuvuefqpvAO519wXA+4ArgZ/1st8bzGy5mS2vq6vLaMyDoREFIiISJqFMBoA7gM3u/vX+Krr7h9x9WfrzXlJJxFVmdnYPde9z91p3rx0/fnymYx6wjj4DdUoGREQkBEKXDJjZImAW8OdD3MXG9PuMjASUBR23CXY11OHuAUcjIiKFLlTJgJldT6qZ//3u3m5m083s0j7qz0lv09UJ6fct2YpzuCqLR1EWL6WprZmDLYeDDkdERApcaJIBM/sA8AXgdmCOmdUClwEXdalzm5mtNbOSdFE18Fkzq06vLwY+BzwPrMhl/INhZtSMmgDAzsN9TaEgIiKSfaFJBoAfANOAJcCy9Ovb3eqUAGWApZdfJNVZ8FEzWwI8Q2o2wmvcPZH1iIdh8qiJAOw4pGRARESCFZqhhe4eH0CdW4BbuizXA5/PZlzZUpNOBnYe6W+whIiISHaFqWWgoKhlQEREwmJYyYCZnZepQApNRzKw87BaBkREJFjDbRn4aUaiKECTRqXmOdjVUEciGeruDSIiMsL122fAzH7S2ypgbGbDKRwlsWKqy8awr3E/dQ37mJQeXSAiIpJrA+lAeCnwYeBIt3IDFmQ8ogIyedQE9jXuZ8fh3UoGREQkMANJBpYAh919afcVZvZixiMqIDWjJrJm9wZ2HN7DOUEHIyIiBavfZKCHRwt3XXdZZsMpLJ0jCjTxkIiIBGjQHQjNbFI2AilER0cUKBkQEZHgDGU0wRMZj6JAdUxJrJYBEREJ0lCSAeu/igzE+LJqYpEY+5sO0tzWHHQ4IiJSoIaSDOiZuxkSiUSYVJGab2CHJh8SEZGAaDrigHX2GziiWwUiIhIMJQMB6+w3oGcUiIhIQIaSDGju3AzSMwpERCRog04G3H1uNgIpVJMrNdeAiIgEKyO3Ccys0szeY2ZnZmJ/haSmS8uAu/pmiohI7g0pGTCzO8yszszONbMyYBnwA+A5M/uzjEY4wlUWV1BRVE5TezMHmg8FHY6IiBSgobYMLATOcPdlwAeBMcA0YAZwU0YiKyAdnQg1E6GIiARhqMlAk7vvTX/+APBdd9/r7ruAxsyEVjj0jAIREQnSQJ5a2JNRZnYSMB14C/AJADOLAmUZiq1gHJ2WWCMKREQk94aaDNwFvEaqZeEH7v6ymZ0P3AmszVBsBUMtAyIiEqQhJQPu/qCZ/S8w0d1Xp4u3ALcC6zMUW8HQ0wtFRCRIQ20ZwN13Aju7LO8AdmQiqEIzqWI8hrH7yF7aE+3EokO+LCIiIoOm6YhDoChWxMSKcSQ9qVsFIiKSc0oGQmLq6BMA2Hxge8CRiIhIoQlNMmBmV5vZo2b2GzN7zsz+x8zeNIDtKs1ssZktM7OVZnanmeVdO/vUqskAbDmoZEBERHJrwMmAmd1lZpeZWVGWYlkM/NDd3+bu5wMvAL8xs4kD2C7q7ucC5wMXA1/JUoxZ09EysOWgul2IiEhuDaZl4Ajwj8BeM/uVmd1oZidmMJal7v5gl+WvA+OAt/e2QfpZCO9Jx4W7t5Ia9rjIzCoyGFvWTa3qSAbUMiAiIrk14GTA3b+YfmLhTOCXpL6k15rZmnTT/AIzG/JtB3e/tltRU/q9uI/NLgWaOXZug2VAKXDRUGMJwqTy8RRF4+xr3M+R1oagwxERkQIylEcY73T377j7e0n9cr8ZiALfBvaZ2U/M7KNmVj3M2OaT+qL/VR91pgO7/djH/e3qsi5vRCIRplTWALBVtwpERCSHhtWB0N3b3P237n6Lu88C5gJPA38C/PlQ92tmBnwJ+KK79zVHbznQ0q2sY/m4aZHN7AYzW25my+vq6oYaXtZ03io4oGRARERyJ6O97t39DeBb6ddw3AFsdvev91OvgeNvI3QsH/fAJHe/D7gPoLa21ruvD9rRToTqNyAiIrkTuiF4ZrYImAW8dwDVNwETzMy63CqYlH7fmIXwsmrq6PTwQs01ICIiORSaeQYAzOx64Erg/e7ebmbTzezSPjZ5klRnwdldympJdT58JnuRZkfnbYJDOzi2G4SIiEj2hCYZMLMPAF8AbgfmmFktcBldRgWY2W1mttbMSgDc/SXgIeAz6fVxUh0a73L3Izk+hWGrKqlkdPEomtqa2dtYH3Q4IiJSIIaUDJjZ6ZkOBPgBMA1YQmp44DJSIxS6KiHVMdC6lH00HdMy4HngD6SenpiXjs5EqE6EIiKSG0NtGVhlZneb2ZhMBeLucXe3Hl7/0KXOLe4+3d2bupQdcvePuPu57n6Ou3/W3dszFVeundj5jIJtAUciIiKFYqjJwJtJ3ad/1cz+2syiGYypoGlaYhERybV+kwEzOyt9L76Tu69x90uB64FPAmvM7B1ZirGgnJTuRLhVIwpERCRHBtIy8Hugx2cQuPvDpFoIvgf8KP3UwWz0JygYUyprMIwdh3fTlmgLOhwRESkAA0kGZgGb+1hfBqwglRBcDrxoZt80s9EZiK/gFMeKmFQxnoQn2X5od9DhiIhIAeg3GXD3re6e6FpmZovM7AEzewXYBzwCnAvcTerWwenAOjM7Lwsxj3gndo4o0K0CERHJvqF2IPw0EAfuJTUPwGh3n+/uf+Pu33f3twPfBO7PUJwF5SR1IhQRkRwa0nTE7t5jH4JuvkvqGQMySEcfWKThhSIikn3ZnIGwDrgki/sfsTS8UEREcilryYCnPJ2t/Y9kE8vHURSNU990gCMtDUGHIyIiI1xonk0gR0UiEU6qmgLApv1bAo5GRERGOiUDITVj7DQAXqt/I9A4RERk5FMyEFKnjD0JgNfq+5riQUREZPgymgxooqHMmZFOBjaqZUBERLIs0y0D+83sFTN7b4b3W3AmjZpAabyE/U0HqW86EHQ4IiIygmU6GXg7cA9wp5m9P8P7LigRi3RpHdCtAhERyZ4hJQO9PYzI3Z9y93tITU28fjiBCZzS0Ylw3xuBxiEiIiPbUFsGVpnZ3WY2pqeV7r7f3V8cRlzC0U6EahkQEZFsGmoy8GZSjy5+1cz+2syiGYxJ0jqTgf2bcfeAoxERkZFqSMmAu69x90tJPaHwk8AaM3tHRiMTqkvHMLqkkobWRnYfqQs6HBERGaGG1YHQ3R8m1ULwPeBHZvZob/0JZPDMrLMToeYbEBGRbMnEaIIyYAWphOBy4EUz+6bmHMiMUzQToYiIZNmQHmFsZotIjRg4FzgFaAVWA3en3z8ErDOza939+UwEWqg0vFBERLJtSMkA8GngWeBe4Dlghbu3dln/fTP7W+B+UrcRZIimp5OB1/dvIZFMEI2or6aIiGTWkJIBdz9xANW+C9wxlP3LUZXFFUwor2ZPwz62HdrZ+TRDERGRTMnmg4rqgEuyuP+CMUOTD4mISBZlLRnwlKcHu52ZFZnZ18ys3cymDaD+ejNb0u1185CCDqmOToTqNyAiItkw1D4DWZH+8v8v4BVgoDfHd7n7wmzFFAaaiVBERLIpm7cJhqIC+DCp/gaSNn3MiZgZWw5upzXRFnQ4IiIywoQqGXD3te7+WtBxhE1JvIQplTUkPMmm+i1BhyMiIiNMqJKBISo3s/vNbGm6v8AXzawk6KAy7fRxpwDwct2rAUciIiIjzUhIBjYA97r7AuB9wJXAz3qqaGY3mNlyM1teV5dfc/2fMf5UQMmAiIhkXt4nA+7+IXdflv68F/gScJWZnd1D3fvcvdbda8ePH5/jSIdnVjoZWL93I4lkIuBoRERkJMn7ZKAHG9PvMwKNIsPGllUxsWI8ze0tvHFgW9DhiIjICJLXyYCZzTGz67sVn5B+H3E97TpaB9bt0a0CERHJnLxKBszsNjNb26WDYDXwWTOrTq8vBj4HPE/qSYojyhnjU40d6jcgIiKZFLZJh4qAJ4CqdNGPzGyHu1+bXi4h9chkSy+/SKqz4KNm1kRqnoLVwF+4+4i7sT5rwmkAvLz3NZKeJGJ5lcuJiEhIhSoZSD/5cGEf628BbumyXA98PvuRhcOE8mrGlY1lb2M9Ww/u0EOLREQkI/TTMs903CpQvwEREckUJQN5ZlbnfAOaqFFERDJDyUCeOWPC0cmH3D3gaEREZCRQMpBnaiomUFVSycGWw+w4vDvocEREZARQMpBnzKxzamL1GxARkUxQMpCHOicfqnsl4EhERGQkUDKQh45OPvSa+g2IiMiwKRnIQ1NG1zCqqJz6pgPsbtgbdDgiIpLnlAzkoYhFOkcVrN29PuBoREQk3ykZyFNnT5oFwKqdLwUciYiI5DslA3nq7JrZAKzZvZ72RHvA0YiISD5TMpCnxpWN5cTRk2lub2H9Xs1GKCIiQ6dkII/NrTkTgJW6VSAiIsOgZCCPnZO+VbBayYCIiAyDkoE8dtq4UyiNl7Dt0E7qGvYFHY6IiOQpJQN5LBaJ8qaJZwAaVSAiIkOnZCDPzU3fKli1c23AkYiISL5SMpDnOoYYrt29gbZEW8DRiIhIPlIykOfGllZxUtUUWhKtvFynIYYiIjJ4SgZGgKO3CtRvQEREBk/JwAigfgMiIjIcSgZGgNOqp1MWL2XH4d3sPlIXdDgiIpJnlAyMANFIlLPSDy5asWNNwNGIiEi+UTIwQrx5ylkAPLt1ZcCRiIhIvlEyMELMq5lDUTTOhr0b2de4P+hwREQkj4QuGTCzIjP7mpm1m9m0AdSvNLPFZrbMzFaa2Z1mFstBqKFSEi/hnJo5gFoHRERkcEKVDKS//J8GaoDoADdbDETd/VzgfOBi4CvZiC/s5k89B4Bnt64IOBIREcknoUoGgArgw8B3B1LZzM4E3gP8I4C7twJ3AYvMrCJLMYbW3JozKY4W8eq+1/XgIhERGbBQNae7+1oAM5sywE0uBZqBrgPslwGlwEXAYxkNMORKYsXMmzyHP2xdwbNbV3LN6ZcFHZJkkbvjDt71s3tqOemd5am6qXp0rPej+6DrMqmK3u04xx6X47dJfTi2Xg/x9n0+fa4+epy+Kw3KIKsPWn/nLNKfqopiKsqKsn6cUCUDQzAd2O3H/sXt6rKu4MyfOi+VDGxZoWSgB4lEkubWBC1tCZpb22lJf25tS9Dalky/J2hrT9LanqStPUlbe4L29iRtidRyeyJJe8JJJFJliYTTnkiSTKbeE0knmXQSyVSdRPpzsst7MukkPf1KL7tDIunpL3Ynmf5y73jv+oWf1HeMSEH42LvP5JqLT8n6cfI9GSgHWrqVdSyXda9sZjcANwBMnTo1u5EFZO6k2ZTEitm4fzO7j9QxsWJ80CFlVHsiyeGGVg51vBpbOdLYypHGNo40pV6NTW00trTT0NRGU0s7jS3tNLe009TSTlt7MuhTyLiIAWapd6xz2QzSH4Gjy5h1Ke+yzjqXoKNuR40udbsuW7cK1nUjjt1HTxW6r+++/XGbH7/HnioNyiCrD1p/5yTSl4rSeE6Ok+/JQANQ3K2sY7mxe2V3vw+4D6C2tnZE/rYqihVRO/lN/H7LMp7dupJ3n3F50CENSDLpHDjSwp76RvYebGLfwWb2Hmii/mAz+w+3cOBIMwcOt3C4cXhPZowYFBfFKCmKUlwUpaQoRlE8QlE8SlE8SnE8SjwWoSiWeo/HI8SjEeKxKLGYEY9GiEUjxGLp92iEWNSIRiNEI0YsGiFiRjRqR5cjRiSSWo7Y0c+W/pwqo8vn1FdeJF0n9UWd+pJPLR/9Mjd904hIBuR7MrAJmGBm1uVWwaT0+8aAYgrc/KnzUsnAlhWhSgYSSWdPfSPb646wo+4IO/Y2sKPuCHv2N7Jnf9OAfrVHDCrLixlVXkRleRGjyuKMKiuioqyI8tIYFaVFlJfEKCuNU14Sp6wkRmlx6lVSHKMoFtEXqIhIN/meDDwJ/Aswm6OdCGuBJuCZoIIK2lmTZlEaL+H1A1vZeXgPNaMm5DyGI42tvLbtAJu2H2TzrsNs3nWIrbsO09rHF/7oiiLGjyljfFUp1ZUljB1dQvXoUsZWFlM1qoSqilQSEI3oy1xEJJPyKhkws9uAdwO17t7s7i+Z2UPAZ4CPmFkcuBm4y92PBBhqoIqicc6dfBZLNz/PM1uW8b7ZV2X1eImk8/qOg6zbtI8Nm/fz6tYD7NzX0GPdcaNLOGFCBZPHVTB5fAWTx5UzsbqMCWPKKC3Oq/8dRURGjFD962tmRcATQFW66EdmtsPdr00vl5DqGNj1p+FHgXvMbBmpiYqeAm7NScAhtmDaeSzd/Dz/u+kPXHvGO4hEMjelhLvzxs5DrFi/hzWv7eXlN+ppamk/pk5RLMLJk0czfcpoTq6p5KSaSqZOqsxZZxgRERm4UCUD6UmDFvax/hbglm5lh4CPZDey/HPmxJlMKK9mT8M+Xtz9MmfXzB7W/ppb2lmxfg/LX97Nyg27qT907CCOSdVlzDq5mlknj+W0qWM4ceIoYtGwzWklIiI9CVUyIJkTsQiXTL+QH635FU9t+v2QkoGmlnaWr9vNMy/uYPn63bS0JjrXja0s5pyZEzn7tPGceUo11aNLMxm+iIjkkJKBEeytJ1/AT9b+mhXbX+RA00GqSkf3u4278/Ib9Tz5/BZ+98L2YxKAmSeN4fwza5h3+gSm1VSqV76IyAihZGAEG1M6mnmT57Bs+wsseeO5PocZNjS18eQfN/P4c5vZtudo38szpo3lwrMmM39ODRPGHDePk4iIjABKBka4S0+5iGXbX+A3G3/PNadfRsSOvY9ft7+JX/1uI48/t7mzE+CYUcW87dypXPbmqUweX3DPexIRKThKBka4sybOYlzZWHY37OWlPa8wZ+LpAOzYe4QfP/kKT6/cRiI90f2bZozj6oumc+6sier8JyJSQJQMjHCRSIRLpqf6Djy18fdMKZ3Gj57cwOPPbSaRdCIGF599Au9ZeAqnnjgm6HBFRCQASgYKwFtPvoCfrv1vntu6imceHUNLU4yIwdvOPZEPXDaTSdXlQYcoIiIBUjJQAF7f3ErkyEQSFbtIjN7KedPn8+Erz+CkSZVBhyYiIiGgZGAE23+omf/45Vp+t3o7kdGTKZ65i+oZO/m7a+YRi+rSi4hIinqJjVBLV23j4//4W363ejvFRVE+cvFCplTWcKj1EL/fsizo8EREJESUDIwwza3t3POT1fzTD1fQ0NTGOadP4FufuYRr33oq7zr97QD8cv0TJL3/xwWLiEhhUFvxCLJl1yHu/MFytuw6TDwW4WPvnsMV55/UOVPghSedy4/W/orth3axcscaak84K+CIRUQkDNQyMEI8vXIbn7prKVt2HWbKhAq+fvMC3jF/2jFTBsciUd4581IAHn75Cdw9qHBFRCRElAzkOXfnx09t4J8fWEFrW4JLak/kG4vewsmTe34OwSXTL6SiqJxX9m1i/d7XchytiIiEkZKBPNaeSPKvP32BH/7PeszgY+86k0/96TmUFvd+96ckVswVp74FSLUOiIiIKBnIU43NbXz1O8/zxPObKYpH+dxH3sw1C04Z0LZXnPpWiqJxVu1cy+YD27IcqYiIhJ2SgTzU2NzGrf/+LCs37GF0RRF3fPwC5s+pGfD2lcUVXDL9QgAeWvdYtsIUEZE8oWQgzzS1tPMP//EcG7bsZ+LYMv75kwuYedLYQe/nmpmXEY/E+MPWFWys35yFSEVEJF8oGcgjzS3tfPk/n+PlN+oZP6aU2z9+4ZCfKzCufCzvOO2tAPzwhV9oZIGISAFTMpAnWtoSfPX+53lp0z6qR5dw+19dyMSxZcPa57vPuJzyojJe2vMKq3a+lKFIRUQk3ygZyAOJpHPn95fx4mt7GTOqmNs/fiE144b/pMGKonLeO+sdADzwwi9IJBPD3qeIiOQfJQN5YPGvX2LZut2MKivi9o9fyAnjKzK278tnvIXx5dVsPbSTp994LmP7FRGR/KFkIOSe+uNmHn56I9GI8bmPnsuJE0dldP/xaJw/nfMuAH689hGa21syun8REQk/JQMh9tKmfXzrZy8A8PH3nsWcU8Zl5TgXTJ3H9DFT2d90kP/e8JusHENERMJLyUBI7a5v5I7Ff6Q94VyzYDqXn39S1o4VsQgfOutaAB56+TF2HanL2rFERCR8QpcMmNk1ZrbMzJaa2TNmVttP/fVmtqTb6+ZcxZsNbe0Jbv/u8xxqaOWcmRP4i6tnZ/2YZ06cyUVTz6U10cZ9yx7QUEMRkQISqkcYm9k84EHgze6+zsyuBh43s9nuvquXzXa5+8KcBZkD33/0ZV7fcYiaceV89sO1RKO5ydk+OvdPeGHXOtbu2cD/vv6HzlkKRURkZAtby8DngMfdfR2Au/8a2A3cFGhUObT6lT08/PRGIhHjlg/Oo7w0nrNjV5aM4qNz3w/AD1b/nP1NB3N2bBERCU7YkoFLgeXdypYBlwUQS84damjlX/5rFQDXvX0mp00dk/MYLjrpXObWzKahrYn7V/4458cXEZHcC00yYGZjgdHAzm6rdgHT+9i03MzuT/cxWGJmXzSzkqwFmiXuzr/+dDX1h5o5Y9pY3ve20wKJw8z42LzrKIkV8/y2Vfxx2+pA4hARkdwJTTIAdEyp132gewvQ17y7G4B73X0B8D7gSuBnPVU0sxvMbLmZLa+rC1eP+af+uIVn1+yktDjG31x3DtGIBRbLuPKxXPemdwPwH8sf1O0CEZERLkzJQEP6vbhbeTHQ2NtG7v4hd1+W/rwX+BJwlZmd3UPd+9y91t1rx48fn5moM2DfwSb+45drAPira9805IcPZdLbZyzgzAkzOdhymLuf/Y6mKhYRGcFCkwy4ez1wAJjUbdUkYOMgdtVRd0YGwsqJ+3/1Ek0tCc6bPYm3zpsSdDhAau6BT57/51SVVLKu7lV++tKvgw5JRESyJDTJQNpTQPd5BWrT5ccxszlmdn234hPS71syHFtWvPBqHUtXb6coHuVj756DWXC3B7qrKh3NzfP/EjPjF+seY9XOtUGHJCIiWRC2ZOBrwOVmdgaAmV0J1ADfSi/fZmZru3QQrAY+a2bV6fXFpIYnPg+syHXwg9XWnuTfH3oRgPdfeuqwH0mcDbMnnMb/OfOdAPzrc4vZ21gfcEQiIpJpoUoG3H0F8EHg+2a2FPgCcHmXCYdKSHUm7Pj5/CKpzoKPmtkS4BlSow+ucffQ3+R+5Hcb2br7CDXjyrl2YXjvarz7jMs5a9IsDrc28C9/+E9aE21BhyQiIhlkhTrtbG1trS9f3n1Kg9zZe6CJj9/5G5pbE/zDx85n3ukTA4tlIA61HOFvn7iDfY37OW/KXD41/3oikVDlkiIi0r8e70XrX/OAfOdXa2luTTB/Tk3oEwGAyuIKPnfxTZTFS3l+2yoWr/qpnl8gIjJCKBkIwPo36vn9Czsoike5/l1nBh3OgE2tOoHPXPRXxCIxHnttCb9a/2TQIYmISAYoGQjADx97GYB3LZjOhDHh6zTYl9kTTuMT530UgAdefIilbzwfbEAiIjJsSgZybM3Gvbzw6l7KSmK8J8SdBvtywdR5fOTs9wFw7x+/zx+2BNf3QkREhk/JQA65Ow88th6Ady84hVFlRQFHNHRXzXwb7znjChKe5O7n7ud/N/0h6JBERGSIlAzk0OpX6nhp0z4qSuNcs+CUoMMZtg/MuYb3n/lO3J17l/2Ax15dEnRIIiIyBEoGcqRrq8C1b51BeWk84IiGz8x43+wr+bOz3wvA/St/zMMvPx5wVCIiMlhKBnJk+cu72bBlP6Mrirj6or6eyJx/rp55KR+bdx2G8eCLD3Pf8gdpT7QHHZaIiAyQkoEccHceeDzVKvC+S06ltDgWcESZd9mMi/nr8/+ceCTGUxt/x1eW3MXB5kNBhyUiIgOgZCAHVm7Yw8ZtBxkzqph3XHBy0OFkzUUnncuXL/k0Y0pHs37vRv7uya/x+v6tQYclIiL9UDKQA79augmAaxacQnE8GnA02TWjehpfu+xznFY9nX2N+/nSb/6Jx199WrMVioiEmJKBLNu6+zArN+yhKB7l8vNPCjqcnBhTOpq/f+siLpl+Ia2JNr6z8kfcsfRfqW86EHRoIiLSAyUDWfbI71KtApfUnpjX8woMVjwa56/O/RB/c8HHqCgq54Vd6/j0Y1/VBEUiIiGkZCCLDje28tsVqXvm77xo5PYV6Mv5J57D16/4EnNrZtPQ2shdz36H/7f0W+w8vCfo0EREJE3JQBY98dxmWloTzD1tPFMnVQYdTmDGlI7m7y6+ievn/Sml8RJW7VzL3zz2FR544SGa25qDDk9EpOApGciSRCLJr595HWBEzDY4XGbG22cs4O4rv8zCk+eTSCb45fonuPl//oGnNv5O8xKIiARIyUCWPLd2F3sPNHHC+HLOmTkh6HBCo6qkkv/75j/j9ks/yyljT2J/00HuW/4gn3z075UUiIgERMlAlvxy6UYA3nnRdCIRCzia8Dm1+mRuv/SzLJr/l0yprGFvYz33LX+Qmx/9ex595bc0tjYFHaKISMGwQh3/XVtb68uXZ6dn+2vbDvCpf3ma8pIY37318hE542AmJZNJntu2kp++9N9sP7QLgJJYMQumnccVpy5kSmVNwBGKiIwYPf461bdUFvzmj1sAuOTcqUoEBiASiXDB1FrOn3IOy3a8wGOvLuGlPa/wxGtLeeK1pZwx/lQuPunNzD/xHMqLyoIOV0RkxNE3VYa1tSd5etV2AN5We2LA0eSXSCTCeVPmct6UuWw5sJ3HXl3C7zb/kZfrXuXlulf57sofM2/ym7hg6jzOmjSL0nhJ0CGLiIwISgYybPnLuzjc2Mq0mkqmnzA66HDy1tSqE7jh3A/yobOv5Y/bVvO7zc+zdvcrPLdtJc9tW0ksEmPOxJnUTj6LcyafSXXZmKBDFhHJW0oGMuw3y1KTDF1SeyJm6jg4XGXxUhaePJ+FJ89nX+N+ntmyjGXbXuCVfa+zaudLrNr5EqyAmlETmDPhdM6cOJPZE05jVHFF0KGLiOQNdSDMoINHWvjIlx/HgcVfejtjKtWMnS0Hmg+xcscalm9/kZf2vEJT+7GTF9WMmsBp1dOZOW46M8aezJTRNcQiI/shUSIiA6AOhNm2dNV2Ekmn9oyJSgSyrKqkkkumX8gl0y+kPZlgU/1m1uxez5rd63m1/g12Ht7DzsN7ePqN5wCIRWKcWFnDtDEnclLVCUyprGFKZQ1jSkerBUdECl7okgEzuwb4EtAERIGb3b3Xn/BmVgl8E5idrv8k8AV3z/nsNb9dnh5FoI6DORWLRDlt3HROGzed986+kvZEO28c2MYr+zbxyr7X2VS/mV1H6nj9wFZeP7D1mG1L4yWcMGoSEyvGMbFiPJMqxjOxYhzjy6oZUzqaqFoTRKQAhCoZMLN5wIPAm919nZldDTxuZrPdfVcvmy0GGtz9XDMrApYAXwE+n4uYO2zedYjXth2kvDTOebMn5fLQ0k0sGmNG9TRmVE/jynRZU1szmw9s4/X9W9lycAfbD+1k+6FdHG5t4LX6N3it/o3j9mNmjC2tYlzpGMaUVjGmdDRjSkdTVVJJVUkllcWjGF0yisriCuLReE7PUUQkk0KVDACfAx5393UA7v5rM9sN3ESqteAYZnYm8B7gTen6rWZ2F7DYzO5w9yO5Cvy36Y6DF599AkVx/ZoMm9J4CaePn8Hp42ccU36o+TDbD+9i95G97DpSx+4jdew5spe9jfs50HyIfY372de4v9/9l8SKqSgqZ1RRORXFZZTFyyiPl1JWVEZZvJSyeAmlsRJK46lXSayYklgxxbFiSqJFFMeKKYrG1RIhIoEIWzJwKXBnt7JlwGX0kAyk6zcDa7vVLwUuAh7LQozHSSSSLFmZHkUwT7cI8kllySgqS0ZxxvhTj1vXnminvukAexvr2d98kP1NR1+HWo5wsOUwh5oPc6jlMM3tLTS3t7C3sX5Y8cQiMYqicYqjRcSjMeLROEXROPFInHg0RiwSIx6JEUt/jkWiR98tSjSSesUiUSIWSb9HiVqEaCTS+TmSfqXKOl6Gpd87ygxLlxtGpMvnVDnHvadaVDrqGGAWSfVYsvQy1u0zdNZOd9+w9IfO3hxd66eXj1l/zLru++jWJ+S4Retrdeexeqt/XPU+1/a20fD6rajXy8iV+jvN/o+E0CQDZjYWGA3s7LZqF/COXjabDuz2Y4dE7OqyLidWv1pH/aEWJo8r5/RpGu8+UsSiMSZUjGNCxbg+67k7Te3NHGlp4EhrA4dbG2hobaKxrZGG1iYa2hppamumqb2Z5raW1Hs6eWjpfG+lNdFGe7Kd9mQ7jW16NoOIwEfn/glXnnZJ1o8TmmQAKE+/t3QrbwF6m4O2vJf69LSNmd0A3AAwderUoUXZgz31jZQURTW3QIEys/StgFIm0Hfi0Bd3py3ZTmuildb2NtqSbbQl2mlNtKUThTbaku20JdrTSUMi/UotJ5JJEp7oXJf0JIlkgkQyQbsnSLqT9CTJdL2kJ0m6pz+n1nlHnfRnx0m64+m6jqfKO+rjeDp292TqMw4OSZKkPnrnvjx1osfU61yTTunTtfCj/2GObptePmZ9l226r3O6DZ0+btH7Wt25v97qH1e9z7W9bTS84d2FOTi8cEQsN88TDFMy0JB+L+5WXgw09rFNT/XpaRt3vw+4D1LzDAwtzOO944KTWTjvRJJJ/VnK0JkZRenbAhQFHY2IFJLQJAPuXm9mB4DuXfEnARt72WwTMMHMrMutgo7te9smK/RAIhERyVe5aX8YuKeA2m5ltenynjxJqrPg7G71m4BnMh6diIjICBS2ZOBrwOVmdgaAmV0J1ADfSi/fZmZrzawEwN1fAh4CPpNeHwduBu7K5bBCERGRfBaqtm13X2FmHwS+b2YdMxBe3mXCoRJSHQO79tL7KHCPmS1L138KuDV3UYuIiOQ3PahIRESkcPQ45C1stwlEREQkx5QMiIiIFDglAyIiIgVOyYCIiEiBUzIgIiJS4JQMiIiIFDglAyIiIgVOyYCIiEiBK9hJh8ysDtic4d2OA/ZmeJ9BGCnnATqXsBop5zJSzgN0LmGUjfPY6+5XdC8s2GQgG8xsubt3f9BS3hkp5wE6l7AaKecyUs4DdC5hlMvz0G0CERGRAqdkQEREpMApGcis+4IOIENGynmAziWsRsq5jJTzAJ1LGOXsPNRnQEREpMCpZUBERKTAxYIOIF+Y2TXAl4AmIArc7O7L+6hfCXwTmJ2u/yTwBXdvz0G4fRrCuawHdnUrfsjd785elANjZkXAV4BbgBnu/kY/9UN5XYZwHqG7JmZ2NfB/gWKgHNgP/K27v9jPdqG7JsM4lzBelwXAImAsqf++VcB/9hWTmcWBO4G3AA6sAha5e0O24+3LEM9lSQ/Fz7r757IQ4pCY2SeAe4C3uvuSPupl72/F3fXq5wXMA44As9LLVwP7gEl9bPML4Afpz0XAH4A78vRclgQddy9xTQOeBb5H6h+saQPYJnTXZYjnEbprQmo89HVdlr8G1AET8/CaDPVcwnhdvg3c2mX5LCABXN3HNt8A/pfUD0YDfgo8mKfnErpr0i2+yaTmvHFgYT91s/a3Evh/iHx4AT8Dft6tbB3w1V7qn5m+sHO6lL0faAQq8ulc0uuXBH0N+vjvPANYOJAv0bBel8GeR1ivCfCLbsvj0+fz4Ty8JoM+lxBfl1nAqG5l+0j90u+p/higFXhnl7I3p89/Rj6dS1ivSbf4fg7c2F8ykO2/FfUZGJhLge7N6MuAy/qo3wys7Va/FLgo49ENzmDPJbTcfa27vzaITUJ5XYZwHqHk7td2K2pKvxf3sVlYr8lQziWU3H2dux8GMLOImX0MaCH1a78nbwHiHPvvxCpSv8AvzWas/RnCuYSamb0TaAMeH0D1rP6tqM9AP8xsLDAa2Nlt1S7gHb1sNh3Y7enUrUv9jnWBGOK5AJSb2f2kfr0mgaeAf3b35qwEmj2hvC5DlA/XZD6pf7x+1UedfLkmAzkXCPF1MbMvAp8A6oEr3X17L1Wnk/oF2tn3wd3bzGwfIbkmgziXjvp3A3NJ3fJ4DrjN3Q9mPdC+YyoHbgcuZ2BJZlb/VtQy0L/y9HtLt/IWoKyPbXqqTx/b5MJQzgVgA3Cvuy8A3gdcSep2Q74J63UZilBfEzMzUp1Uv+jue/qoGvprMohzgRBfF3e/Dagh9QX0tJld0EvVcqCt25cO9P/vRM4M4lwAXgD+J31NriCVqC1Jd5IM0leBb7t79x9nvcnq34qSgf519J7tnrkVk7pX09s2PdWnj21yYSjngrt/yN2XpT/vJfUP41VmdnY2gsyisF6XQcuDa3IHsNndv95PvXy4JgM9l9BfF095AFhKqlNkTxqAeDoJ6qrPfydybYDngrvf7O6PpT83kBqxczZwVS7i7ImZnQOcR6pD5EBl9W9FyUA/3L0eOABM6rZqErCxl802ARO6/TF1bN/bNlk3xHPpSUfdGRkIK5dCeV0yJDTXxMwWkero9ecDqB7qazLIc+lJ4NclPWy1u3Wkhqf1ZBOp5vSJXfYRA6oJ+JoM4Vx68jrpzpAZCWporiJ1r/+36aGPP0qX32VmS8xsZg/bZPVvRcnAwDwFdH9yVG26vCdPkrrQXf8HrSXVCemZjEc3OIM6FzObY2bXdys+If2+JcOxZVuYr8uAhfmapOO6Eni/u7eb2XQz66vTWWivyWDPJcTXZUUPv/InA73dZ3+a1GiCrv9OzCU1rr23f/NyZVDnYmYTzOwL3YprSCU7gV0Td/+qu5/j7gvdfSHwgfSqRemyDT1slt2/laCHVeTDi9TY/MPAGenlK0l1XJmUXr6NVA/Pki7b/AL4XvpzHPg94ZlnYMDnQmq42ytAdXq5GPg1qU440aDPp0uMxw3Jy6frMpjzCOs1IfUP2uukeqPXpl83Av+Qb9dkKOcS4uvyBvDXXZbnkfoCuTm9fCOpce41Xep8A/gNR+cZ+DHhmGdgUOdCag6PfR1/U6R+AP97+toGOsy723lNo9vQwlz/rWg0wQC4+woz+yDwfTPrmLXvcnfv6MlZQqoDR9eM9aPAPWa2jKMZ9a25i7pnQziXF0l1gHo0Xb8CWA38hbsnchp8N+kmwydIzUIG8CMz2+FHh4XlxXUZwnmE9Zr8gNSXx5Ju5V/u8jkvrglDO5ewXpfPA9eb2XWkhgeWAp8G7k2vL0qXRbts83ekZiB8ntSX1Grg5hzF25fBnssu4OvAf5lZM6lOeJuBt7n7kVwG3hszuws4P714l5m96u5/Qo7/VvSgIhERkQKnPgMiIiIFTsmAiIhIgVMyICIiUuCUDIiIiBQ4JQMiIiIFTsmAiIhIgVMyICIiUuCUDIiIiBQ4JQMiklNmdruZfTPoOETkKCUDIpIRZvZPZvb4AKpeBTya7XhEZOCUDIhIprwZ+GNfFcxsMqlHxy7JRUAiMjB6NoGIDEv6QUtHSD1FrcPL7j6rh7rXA+9296v72N91wKfS+ywBlrj75zIbtYh0pWRARIbFzCKknne/HDiP1HPiW9x9fw91fwE86e73dl+XXj85vf1p7r7JzMYD6929OmsnICK6TSAiw+PuSaAGOAwsc/ddvSQCRcCl9N1fYCKpR7NOS++7Drgy0zGLyLGUDIhIJswFXvC+mxovBra6++Y+6qwGfgA8ZWa/NbMbgLWZC1NEeqJkQEQy4WxgVT91rqSfUQSe8mfAHGAFcDuw2syqMhCjiPRCyYCIZMJZwIv91LkS+O++KpjZCWY2391fcvfPALOBycDbMhOmiPREyYCIZEIMON3MJvf0K97MTib1pf77fvZzKnCnmcXSyxHAgFczGKuIdKNkQEQy4QvAB4BtwP/rYf1VpEYRtPezn/XARuBZM1sCPALc5O79tTqIyDBoaKGIZJ2ZPQr8zN3vDzoWETmeWgZEJBeWkPqVLyIhpJYBERGRAqeWARERkQKnZEBERKTAKRkQEREpcEoGRERECpySARERkQKnZEBERKTAKRkQEREpcP8fQSBzLppFoEoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(time,y_analytical(time,0),label='$y$')\n",
    "plt.plot(time,dy_dt(y_analytical(time,0)),label='$\\\\dot{y}$')\n",
    "plt.legend()\n",
    "plt.xlabel('$t$ / s')\n",
    "plt.ylabel('$y$, $\\\\dot{y}$ / s$^{-1}$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to get the same result with a forward Euler scheme. The forward Euler scheme proceeds fully explicitly with the update formula\n",
    "\n",
    "$y_{n+1} = y_n + f(y_n) (t_{n+1} - t_n)$\n",
    "\n",
    "We define the forward Euler update as\n",
    "\n",
    "$\\Delta y_\\text{FE} = y_{n+1} - y_n = f(y_n) \\Delta t$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Forward_Euler_Update(y_n,dt):\n",
    "    return dy_dt(y_n)*dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now need to define the time stepping itself. So we choose the same starting value (initial condition) of $y_0 = y(t=0) = 0$ to initialize our solution field and then select a time step size $\\Delta t$.\n",
    "Now, for every time step we append to this array the value of the next time step by applying the forward Euler update rule. Each new time step value then serves as the initial condition for the following time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrate(dt):\n",
    "    #Startwerte\n",
    "    y = np.array([0.])\n",
    "    times = np.array([0.])\n",
    "    t_end = 4.\n",
    "    while times[-1]+dt < t_end: #repeat the loop as long as the final time step is below the end point\n",
    "        times = np.append(times,times[-1]+dt) #here define the next time point as the previous time point plus the time increment dt\n",
    "        y = np.append(y,y[-1]+Forward_Euler_Update(y[-1],dt)) #here we calculate the solution at the new time step\n",
    "    return times, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "db57f62a6389438d87e4737951e8396b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(BoundedFloatText(value=0.1, description='$\\\\Delta t$ / s', max=1.0, min=0.001, step=0.1)…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ipywidgets import widgets\n",
    "from ipywidgets import interact\n",
    "\n",
    "@interact(dt=widgets.BoundedFloatText(value=0.1,min=1.e-3,max=1,step=0.1,description='$\\Delta t$ / s'))\n",
    "\n",
    "def plot_FE(dt=0.1):\n",
    "    time_n, y_n = integrate(dt)\n",
    "    fig, ax = plt.subplots(ncols=2,figsize=(18,6))\n",
    "    #solution\n",
    "    ax[0].plot(time,y_analytical(time,0),label='analytical')\n",
    "    ax[0].plot(time_n,y_n,label='Forward Euler',marker='d',ls=':')\n",
    "    ax[0].legend()\n",
    "    ax[0].set_xlabel('$t$ / s')\n",
    "    ax[0].set_ylabel('$y$')\n",
    "    #error\n",
    "    ax[1].plot(time_n,y_n-y_analytical(time_n,0),ls=':',marker='d',label='absolute')\n",
    "    ax[1].plot(time_n[1:],(y_n[1:]-y_analytical(time_n[1:],0))/y_analytical(time_n[1:],0),\n",
    "               ls=':',marker='d',label='relative')\n",
    "    ax[1].legend()\n",
    "    ax[1].set_xlabel('$t$ / s')\n",
    "    ax[1].set_ylabel('error')\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tasks:\n",
    "\n",
    "- How does the peak relative error change when you divide the time step size by 10?\n",
    "- How does the solution change when you increase the rate constant by a factor of 2 and why?\n",
    "- How does the solution change again for $k=4$ and time step sices of $\\Delta t= 0.4$ and  $\\Delta t = 0.5$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the stability analysis, the substitution $z = y_\\text{eq} - y$ and $\\bar{k} = -k$ yields our standard test problem $\\dot{z} = \\bar{k} z$ for which we know the stability criterion $|1 + \\bar{k} \\Delta t| < 1$. The two conditions produce $\\bar{k} < 0$ (inaccordance with our decay problem) and $\\Delta t < 2/|\\bar{k}|$."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit 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
}