02_Backward_Euler_Picard.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 2 - Backward Euler, Picard iterations\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 = 2."
   ]
  },
  {
   "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": "iVBORw0KGgoAAAANSUhEUgAAAgwAAAF9CAYAAABoEDD2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABQWElEQVR4nO3dd3xc1Z338c9vZtSbZVm2ZFvuBRdwQQab6oApAUJ62VSyS2DzkATyBJKQEHZJSBbSFpInIWE3iUMIARIgCQRDaKbYBtw77r3KkiWrl5nz/DEjWR5kq3hGc0f6vnnNa+aee865v+uLND/de+655pxDRERE5FR8iQ5AREREvE8Jg4iIiHRKCYOIiIh0SgmDiIiIdEoJg4iIiHRKCYOIiIh0KpDoALzsyiuvdM8991yiwxAREektdrIVOsNwCkeOHEl0CCIiIp6ghEFEREQ6pYRBREREOqWEQURERDqlhEFEREQ6pYRBREREOqXbKkVEpN8KhUIcOXKEyspKgsFgosOJu/T0dIYPH05KSkq32yphEBGRfmvv3r2YGaNGjSIlJQWzk05DkPScc5SXl7N3715Gjx7d7fa6JCEiIv1WbW0tw4YNIzU1tU8nCwBmRkFBAQ0NDT1q75mEwcyuMbNnzewlM3vTzBaY2VldaJdrZvPNbKmZrTCze80sEFWn2Mz+ZmZLInVujd+eiIhIMvH5PPNVGHenkxR56V9pPvCwc+5S59xsYDXwkpkN6UI7v3NuFjAbuBD4butKM/MBTwNrnHNzgPcAXzSzG2K/CyIiIn2TlxKG15xzj7Rb/gkwCLj8ZA3MbCrwQeCHAM65JuA+4BYzy45UuwqYHukP51wV8GvgDuvr559ERERixDMJg3PuQ1FF9ZH3tFM0mwc0AOvalS0FMoAL2tXZ5pyrjKpTAkzsabzd9eaeFfxm+aPsPXagtzYpIiISM55JGDowh3Ay8PdT1BkDHHLOuXZlB9uta32P/paOrtPGzG4ws2VmtqysrKz7UZ/EW3tX8vzWV9l8ZEfM+hQREektnkwYIpcKvgPc4Zw7fIqqWUBjVFnrcmY36rRxzj3onCt1zpUWFhZ2L/BTGJE3DIDdlXtj1qeIiEhv8WTCAPwA2OWc+0kn9Wp59yWL1uW6btSJu5EDIglD1f7e2qSIiPRhw4cP56c//ekJZWvXriU9PZ0NGzbEfHuem7jJzG4BJgMf7kL17cBgM7N2lyWKIu/b2tW5MqpddJ24GxFJGHZV7sU51+fv9xURSVbv+9rfErLdp3/y/m7VnzNnDkuXLj2h7JZbbuH6669n8uTJsQwN8NgZBjO7nvBdDR9zzrWY2Rgzm3eKJi8QHuA4pV1ZKeEBk4va1RlnZgOi6uxxzm2KWfCdKMjIJzMlg+qmWiobjvXWZkVEpI+KThj++te/snLlSu666664bM8zZxjM7BPAt4HrgDMjf4GfDRQDL0bq3A18ACh1zjU459ab2VPAbcDnzCwFuBm4zzlXE+l6AbAK+CrwH2aWC9wA3N07exZmZowcMIyNZVvZXbWP/Iy83ty8iIh0UXf/0k+U2bNn87WvfY2KigqysrK49dZbufPOOykoKIjL9jyTMAB/IBzPwqjy9qlSOuGBiu3P518H/NzMlgJ+wsnFna0rnXMhM7sW+JWZLYn08aBz7sFY70BnSvKGsrFsK7sq9zGtKPani0REpP84++yzSU1NZdmyZaxcuZJAIMBNN90Ut+15JmFwznX66Czn3K3ArVFlx4DPddJuP3DtaQUYAyPzhgOwu2pfgiMREZFkl5aWxowZM3j66af5/e9/zyOPPNKjp1B2lafGMPR1bXdKVCphEBGR0zdnzhx+8YtfMGfOHK655pq4bksJQy8qyRsKwN5jB2kJ9f3nrouISHxNnz4dn8/3rtsr40EJQy/KSElncFYBLaEWDlafaj4qERGRzj388MPceOONTJkypfPKp8kzYxj6ixEDhnO4tpxdVXsZnlec6HBERCTJhEIhysrKmD9/PuvWrePxxx/vle3qDEMvGxG5LLG7UjM+iohI97322msUFxczf/58nnjiCfLz83tluzrD0MtaBz7u0p0SIiLSA3PnziUUCvX6dnWGoZe1PoRqj+6UEBGRJKKEoZcVZReS4k+hrK6Cuqb6RIcjIiLSJUoYepnf56ckNzzYURM4iYhIslDCkACtlyWUMIiISLJQwpAAxx91rYRBRESSgxKGBGi7tbJKt1aKiEhyUMKQAO2fKeGcS3A0IiIinVPCkAB56bnkpedS39JAWV1FosMRERHplBKGBDk+46PGMYiIiPcpYUiQkXmtAx/3JjgSERFJZrt37yY7O5uVK1fGdTuaGjpBRuWXALCjck+CIxERkWQ2dOhQVq1aRUlJSVy3o4QhQUa3JgxHlTCIiEjPBQIBxo0bF/ft6JJEggzNGRKeIrq2nJqm2kSHIyIickpKGBLE7/MzKjKOYedRjWMQERFvU8KQQKN0WUJERHpo+PDh/PSnPz2hbO3ataSnp7Nhw4aYb08JQwKNzh8BaOCjiIh035w5c1i6dOkJZbfccgvXX389kydPjvn2NOgxgY4PfNyd4EhERKTVxx77YkK2+/jHH+hW/Tlz5vDLX/6ybfmvf/0rK1eu5PHHH491aIDOMCRUSd5QfOZjf/UhGloaEx2OiIgkkdmzZ7Nt2zYqKipobGzk1ltv5c4776SgoCAu29MZhgRK9adQklvMrqp97K7cx4RBYxIdkohIv9fdv/QT5eyzzyY1NZVly5axcuVKAoEAN910U9y256kzDGaWamb3mFmLmY3qpO51ZvaOmS2MejWa2Qfb1avsoM7H474zXaSBjyIi0hNpaWnMmDGDp59+mu9///v8+Mc/JiUlJW7b88wZhkiC8CdgM+DvYrN7nHPz2/UxBlgOPNeuzirn3NzYRBl7o/NLeHXnmxrHICIi3TZnzhzuv/9+LrvsMq655pq4bsszCQOQDXwGGA58tgv1/9pB2eeAR51z9TGMK65Ga4poERHpoenTp+Pz+d51e2U8eCZhcM6tAzCz4V2sX9l+2cyMcKLxsZgHF0ejBoQThj1VB2gJthDwe+aQiIiIxz388MPceOONTJkyJe7b6kvfTu8BapxzS6PKi8zsMaAYaCJ8ZuIB51ywl+PrUEZKOsXZgzlQc5i9xw60jWkQERHpSCgUoqysjPnz57Nu3bq43UYZzVODHk/TdcBvOijfCtzhnLsIuB74CvCzk3ViZjeY2TIzW1ZWVhaXQKNp4KOIiHTVa6+9RnFxMfPnz+eJJ54gPz+/V7bbJxIGM8sB3gc8HL3OOXeNc25L5PNO4F7gi2Y2sKO+nHMPOudKnXOlhYWFcYz6OD25UkREumru3LmEQiE2btzIeeed12vb7RMJA+FxCy855450oe42wADPTHqggY8iIuJ1fSVhuI4OLkeY2aVmdm1U8bDIu2fuYxwdGfi4s3IvoVAowdGIiIi8W9IkDGZ2t5mtM7P0qPJxwCjg+Q6alQC3mVlmpG4e8FXgCefc4TiH3GW56TkUZOTT2NLIwRrPhCUiItLGM3dJmFkq8E9gQKToUTPb75z7UGQ5HcgkfDmhveuAh5xzHf1p/hIwE3jFzBqAHOBl4K7YRn/6RueXUF5/lO1H9zA0tyjR4YiIiJzAMwmDc64JmHuK9bcCt3ZQfscp2uwhfFeE543KL2HZ/jXsOLqbC0bOSnQ4IiL9hnOO8FQ+fZ9zrsdtk+aSRF83duBIALZrimgRkV6TkpJCfX3STA582pqbmwkEenauQAmDR7QmDNsqdmngo4hILxk8eDD79u2jrq7utP76TgahUIhDhw6Rl5fXo/aeuSTR3w1Iz2VQ5kCO1FWwr/ogJXlDEx2SiEifl5ubC8D+/ftpbm5OcDTxl5WVxaBBg3rUVgmDh4wbOIojdRVsq9ilhEFEpJfk5ua2JQ5ycrok4SGtlyW2lu9MbCAiIiJRlDB4yLiCUUB4HIOIiIiXKGHwkDH5IzCMnVV7aQ72/WtpIiKSPJQweEhGSjrDcosIhoLsqtyX6HBERETaKGHwmLZxDBU7ExuIiIhIO0oYPGbcwFGAEgYREfEWJQwe0zbwsVwDH0VExDuUMHjMyLxhBHwB9lUfpK6p/0xXKiIi3qaEwWMC/gCjBgwHYPtRnWUQERFvUMLgQccHPiphEBERb1DC4EEa+CgiIl6jhMGDxhZEnlypgY8iIuIRShg8aGjOEDIC6ZTXH+VofVWiwxEREVHC4EU+87WNY9imyxIiIuIBShg8SjM+ioiIlyhh8KjWCZy2lO9IbCAiIiIoYfCsCQVjANhSvpNQKJTgaEREpL9TwuBR+Rl5DM4qoKGlkd1V+xMdjoiI9HNKGDxswqCxAGwu35bgSEREpL9TwuBhEyOXJTYd2Z7gSEREpL9TwuBhEweFE4bNShhERCTBPJUwmFmqmd1jZi1mNqoL9d8xs4VRr5uj6uSa2XwzW2pmK8zsXjMLxG0nYqgkbyjpgTQO1R6hUhM4iYhIAnkmYYgkCK8CxYC/i80OOufmRr3uj6ozH/A752YBs4ELge/GKOy48vv8jI/cXrmpXGcZREQkcTyTMADZwGeA38WqQzObCnwQ+CGAc64JuA+4xcyyY7WdeJpQEB74qHEMIiKSSJ5JGJxz65xzW2Pc7TygAVjXrmwpkAFcEONtxcXE1jsllDCIiEgCJcW1/FPIMrPfAuOAEPAi8GPnXENk/RjgkHPOtWtzsN06zxtfMArD2H50N03BZlL9KYkOSURE+iHPnGHooU3AA865i4CPAFcBf2m3PgtojGrTupzZUYdmdoOZLTOzZWVlZbGOt9uyUjMZnldMS6iFHUd3JzocERHpp5I6YXDOfdo5tzTy+QjwHeBqM5seqVILpEU1a12uO0mfDzrnSp1zpYWFhXGIuvs0H4OIiCRaUicMHWidEnFc5H07MNjMrF2doqi6njdhUGvCkDQhi4hIH5O0CYOZnWlm10cVD4u8t567f4HwAMcp7eqUAvXAovhGGDvtBz6eOBxDRESkdyRNwmBmd5vZOjNLjxQVAF83s4LI+jTgduAtYDmAc2498BRwW6ROCnAzcJ9zrqaXd6HHirILyUnLpqqxmkO1RxIdjoiI9EOeSRgiszwuJDxPAsCjZvZkuyrphAcqtl5eWEN4gOOzkXaLCN8Bca1zLtiu3XWR/pcSTiYWA3fGZSfixMzaxjHo9koREUkEz9xWGZlUae4p1t8K3NpuuQL4Vhf6PQZ8LgYhJtSEQWNYtn8Nm45s46JR5yY6HBER6Wc8c4ZBTq31QVTvlMV6bisREZHOKWFIEmMHjiLFF2DPsQMca0ya4RciItJHKGFIEqn+FMYXjAZ0lkFERHqfEoYkMnnweAA2HN6c4EhERKS/UcKQRCYXRhKGsi0JjkRERPobJQxJZHzBGPw+P7sq91HTVJvocEREpB9RwpBE0gKpjBs4CofjnTJNEy0iIr1HCUOS0WUJERFJBCUMSaZ14OPGw0oYRESk9yhhSDITC8bgMx/bK3dT11yf6HBERKSfUMKQZNJT0hmbPwLnnB53LSIivUYJQxKaNHgCABt0WUJERHqJEoYkpIGPIiLS25QwJKEzCsdiZmyv2EVDS2OiwxERkX5ACUMSykzJYPSAEoIuxOYj2xMdjoiI9ANKGJLU5NZxDGV6roSIiMSfEoYk1TaOQQMfRUSkFyhhSFKTCsdhZmyp2ElDc0OiwxERkT5OCUOSykrNZGz+SIKhIBvKtiY6HBER6eOUMCSxs4rOAGDtoXcSHImIiPR1ShiS2JlDJgFKGEREJP6UMCSxCQWjSfOnsrtqH5X1VYkOR0RE+jAlDEksxZ/CpMJxAKw9tCnB0YiISF+mhCHJ6bKEiIj0BiUMSe7MIccHPjrnEhyNiIj0VZ5KGMws1czuMbMWMxvVSd0hZvYDM3vDzBaa2Uozu93MAlH1KiPr278+Htcd6UUjBgwlLy2H8vqj7K8+lOhwRESkjwp0XqV3RBKEPwGbAX8XmnwUuAKY65yrNrMSYAWQAdzZrt4q59zc2EbrHT7zMXXIRBbtXsbaQ+8wLLco0SGJiEgf5KUzDNnAZ4DfdbH+YeBHzrlqAOfcHuBx4FPxCc+7WscxrNE4BhERiRPPnGFwzq0DMLPhXaz/eAfF9UBaLONKBmdFxjGsP7yJYCiI39eVEzQiIiJd56UzDLEwh/BZhvaKzOwxM3vNzF40sy+ZWZ/6Rh2UNZDinMHUNzewrWJXosMREZE+qM8kDGY2DxgBfDdq1VbgDufcRcD1wFeAn52inxvMbJmZLSsrK4tbvLHWereELkuIiEg89ImEwcxGAg8A1zrnKtuvc85d45zbEvm8E7gX+KKZDeyoL+fcg865UudcaWFhYXwDj6GzNB+DiIjEUdInDGY2BPgbcL1zbmUXmmwDDBgT18B62ZTBEzAzNh/ZRl1zfaLDERGRPiapEwYzyweeAW53zr0aKbuh3fpLzezaqGbDIu+7eyfK3pGVmsnEgjEEXUhnGUREJOaSJmEws7vNbJ2ZpUeWs4EFwItAmZmVmlkpcGO7ZiXAbWaWGWmTB3wVeMI5d7h39yD+ZhRPBWDl/nUJjkRERPoaz9xWaWapwD+BAZGiR81sv3PuQ5HldCCT8OUEgJuBcyOvb56k25eAmcArZtYA5AAvA3fFfAc8YEbxVP609m+sPLAe5xxm1nkjERGRLvBMwuCcawLmnmL9rcCt7Za/D3y/kz73EL4rol8YOWAY+Rl5HK2vYlflXkbllyQ6JBER6SOS5pKEdM7M2i5LrDigyxIiIhI7p5UwmNm5sQpEYmNm6ziGA+sTHImIiPQlp3uG4c8xiUJiZuqQifh9fjaXb6emsTbR4YiISB/R6RgGM+vomQ0QHnzY4eRHkjiZKRlMGjSOdYc3sfrQBs4fMSvRIYmISB/QlUGP8wg/RbImqtyAi2IekZy2GcVTWXd4Eyv3r1fCICIiMdGVhGEhUO2cey16hZmtiXlEctpmDJ3CH1Y/wcqD6wm5ED7T2FYRETk9nX6TOOc+1FGyEFl3WexDktM1LKeIwqwCqhtr2F7Rpya0FBGRBOn2n55mVhSPQCR2wrdXTgF0e6WIiMRGT85V/zPmUUjMHb+9UgmDiIicvp4kDJpvOAlMGTyRFF+AbRW7OFpflehwREQkyfUkYXAxj0JiLi2QyplDzgBg+X6NTRURkdOj4fN92DnDpwOwdN/qxAYiIiJJTwlDH3b20DMxjLWHNlHXXJ/ocEREJIn1JGEIxjwKiYu89FwmDhpDS6iFVXq2hIiInIZuJwzOuRnxCETio/WyxNu6LCEiIqchJpckzCzXzD5oZlNj0Z/ETumwaQCs3L+O5mBzgqMREZFk1aOEwcx+YGZlZjbLzDKBpcAfgDfN7LMxjVBOS1F2ISPyhlHf0sD6w1sSHY6IiCSpnp5hmAtMcs4tBT4F5AOjgHHATTGJTGLmnOHhswxv71uV2EBERCRp9TRhqHfOHYl8/gTwO+fcEefcQaAuNqFJrMwaNh2AZftWE3KhxAYjIiJJqacJQ46ZjTSz9wAXA/MBzMwPZMYoNomRUQOGU5g5kMqGY2wt35nocEREJAn1NGG4D9gKvAg87JzbaGazgZcBPbzAY8yMWZHBj5rESUREeqJHCYNz7hFgBHC2c+66SPFu4E7gW7EJTWJplmZ9FBGR0xDoaUPn3AHgQLvl/cD+WAQlsXfGoLHkpGaxv/oQe6r2U5I3NNEhiYhIEtHU0P2E3+dvuyyxZM/yBEcjIiLJRglDP3LeiFIAFu9ejnN66KiIiHSd5xIGM0s1s3vMrMXMRnWhfq6ZzTezpWa2wszuNbNAVJ1iM/ubmS2J1Lk1bjvgYVMGTyA3LZv91YfYVbk30eGIiEgS6XLCYGb3mdllZpYar2AiCcKrQDHg72Kz+YDfOTcLmA1cCHy3XZ8+4GlgjXNuDvAe4ItmdkPsIk8Ofp+fc4eHHwWyWJclRESkG7pzhqEG+CFwxMz+bmY3mllJjOPJBj4D/K4rlSPPrvhgJC6cc02Eb/m8xcyyI9WuAqYDP4nUqQJ+DdxhZhbD2JPC8csSy3RZQkREuqzLd0k45+4g/CVbTPhL+Crgh2a2G3gW+AfwhnM9n0rQObcOwMyGd7HJPKCBE+d+WApkABcAz0XqbHPOVUbVKQEmAu/0NN5kNGnQOAak53K4tpxtFbsYVzAq0SGJSB/mnMM5cO0/O/eu5XDdcD3a1kevC1dobdtW3roO3rX+hLYdrG//d1NbWVvBCW/v+iOrdfHk7d69jeh/mxOWo9t30M7x7sIxQ/NIT+vxTY9d1u0tRG6n/A3wGzNLIXwJ4CrgV0Cxmb1AOIF42jlXHstgOzAGOORO/Fc/2G5d6/sBTtS+Tr9KGHw+H7NLZvLcloUs3rNcCYNIN4RCjmAoRHNL+NUSDNESdASDIZqDIVpaQgRDjuaWEMHQ8XUtQde2HGotDzlCwRAtIUcwsj7cf+QVDPcVco5Qa33Xvi7hda11QsfrhEIO59y76ri2z5zwubW+cxyv1/rF7trXBzhxXbjNyZMBib+ffW0uo4fmxX07p5WSOOeaCc/u+DJwa2QMwtXAR4FBwI9PN8BOZAGNUWWty5ndqNMmMrbhBoARI0bEJkqPOa+klOe2LGTJnuV8etoH8Znnxr6KdKolGKKhKUhjUwsNTUEaGsPvjc1BmpqDNEZ9bmqOLLeEaGoO0hz13tQSorklSFNzJCEIhmiOrG8JhiJJgL4Be8osPOuszwCsbdkMLLK+rTzS4F3lBka4wsnatG2r/TLt2kYKjrc/3gdRZe2XW4ujr2QfX39ixeh27WPpqH27SE7sh3eLjiEttatD/k5PTM9hOOd2Ar+IvHpDLZAWVda6XNeuTm4nddo45x4EHgQoLS3tk78dJgwaTUFGPuV1R9lSvoOJg8YmOiTpJ5xz1De2UFvfQm1DM7X1zdQ2NFNX30xtQwt1Dc3UN7ZQF/lc19BCQ1ML9Y2tr+OJQUswMQ9SC/iNlICPgD/88vt9kWU7vuz34Y8sB/w+/D4LL/t8+CLvfr9FyiPrI599ZgT8hi9S5vO1fvYdX7Zwfz47vt5ndsJ6ny/yBd263L4sUm5Gu/qR5ciXUWv58S/6E7/kT1gGzBeVDNAuIeh/w8X6pPhf9Iiv7cBgM7N2lyWKIu/b2tW5MqpddJ1+xWc+5pTM5JnNL7F493IlDNIjjc1BjtU0UVXbyLHapsirkeraZqrrmqiubaK6roma+mZq6pqpqW+itr6ZWP2R7vMZ6al+0lP9pKUGwu8pftJTA6RFPqem+ElLDb+npvjaylIDvsh7uDwlxU9KwEdqwEdKIPy59RXwt372E/Cbvvyk30r2hOEF4L+BKRwf+FgK1AOL2tX5kpkNaDfwsRTY45zb1Iuxesp5I0p5ZvNLvLlnBZ+b/hF8Pl2WEAgGQ1TWNFJxrIGKqgaOVjeGX8caOFrdQFVNE5XVjVTWNFLf2NKjbaSn+snKSAm/0sPvmekBMtNTyEoPkJEeIDMthYy08OeMtACZaeH39LRwYpCRFiAl4NOXt0gvSqqEwczuBj4AlDrnGpxz683sKeA24HORQZg3A/c552oizRYAq4CvAv9hZrmExyjc3dvxe8nYgSMZnFXA4dpyNpRtYeqQiYkOSeIsFHIcrW6grLKesop6yirrOFLVwJHKeo5U1lNeVU9ldWOXzwD4fUZediq5WWnkZqW2vXJaP2emkp2ZSnZmCjmZqWRHkoSAX8mpSDLqUcJgZmc452J+d0FkUqh/AgMiRY+a2X7n3Iciy+mEByq2/7PiOuDnZraU8GRPLxJ+aiYAzrmQmV0L/MrMlkT6eDAyVqHfMjMuGDmLJzc8x2u73lLC0EfUN7ZwsLyW/WW1HCiv5VBFHYfKazlYUUfZ0fpOr/ubwYCcNAbmpjMwN538nDTyI+8DctLIz0knLzuVAdlpZGWk6C98kX7EejJ5j5nVEx4Y+J/OuaMxj8ojSktL3bJlyxIdRtzsP3aQWxbcRXogjQfffy/pgejxo+JFwZDjUEUtew/XsO9wTfi9rIb9ZTUcrY6+IehEuVmpDM7PoDA/k8L8DAoHZDBoQAaD8jIoyMsgPzdNZwBE+reT/hXQ00sS5xAeO7DFzO4CfumcC/awL0mQoblFjC8YzZbyHSzdu5oLR52T6JCkHeccRyob2Hmgip0HjrH7YDW7D1az93A1TS0dnykI+H0UD8qkuCCbooJMigqyGFKQSdHATAbnZ/bK5C4i0jd1+tvDzKYBGyJzLgDgnFsLzDOzDwA/Ivxshq855xbELVKJi4tHncuW8h28tutNJQwJFAyG2FtWw7a9VWzbV8m2veEkoba+ucP6BXnplAzOYfjgbIYNzmZYYfhVMCADv0+XCUQk9rry58YbwDTCtyeewDn3VzN7lvCAwkfNbBHwf+MxvkHi47ySUuav/AtrDr1DRV0lAzMHJDqkPs85x6GKOrbsrmTznqNs3n2UbfuqaGx690m6nMxURg/NZdTQXEYW5TKiKIeSwTlkZaQkIHIR6c+6kjBMBvafYn0msBz4PXATsMbMfgV8J/KgJ/Gw7LQszh56Jm/tXcnru97m/ZMuT3RIfU5LMMS2vZVs3FnBhh0VbNxZQWUHYw2GDMxk7PA8xg4bwJhheYwZlkd+TpoGFoqIJ3SaMDjn9kSXmdktwKzIayzQRPjWxfsj758GNpjZh5xzb8UuXImHi0edy1t7V/Lqzje59ozL9AV1mlqCIbbsrmTttiOs23aEjTsraIg6e5CTmcrEkflMGJHPhBEDGF+ST25W3J4cLyJy2no6AuprwBLgAeBNYHnk0dKtHjKzbwC/JTypknjY9KIp5KRls/fYAXYc3c2YgSMTHVJScc6x93ANKzcfZtXmMtZtO0J944kJwrDCbCaPHsjk0QOZNLqAoYOylJiJSFLpUcLgnCvpQrXfAT/oSf/SuwL+ABeMmMWCLa/w6s63lDB0QUNTC2u2HmHZxkMs23iIsqP1J6wfVpjNWeMGcebYQUwdW0B+bnqCIhURiY143mNVBlwSx/4lhi4edS4LtrzCG7uX8pnpHybg652nnyWTqppG3l5/kMVrD7BmS9kJtzbmZacybXwhMyYUMm38YArzMxIYqYhI7MUtYYg8DOrVePUvsTU6fwQlucXsOXaAlQfWMWvYtESH5AlHjzXwxur9LF67nw3by0+YNnlcyQBmTRpC6aQhjBs+AJ9uZxSRPkyzuAgQnip67ujz+MPqJ3hp2xv9OmE4VtvEkrX7eW3lPtZtO9KWJAT8xvRxhcw5s5hzpxTpMoOI9CtKGKTNxaNn86e1f2PlwfUcqa1gUNbARIfUa1qCIZZvPMRLy/awdMNBWoLhLCHg9zHrjMFcMG0osyYXaf4DEem3lDBIm9y0bM4ZPp3Fu5fx8o5FfGzq+xIdUtztOVTNc2/u5NUVe6mqCd/o4zOYMaGQi2YMZ/aZxWQrSRARiW3CYGZ5mqwpuc0bc0E4Ydi+mA9Pvgp/Hxz82NwSZPGaAyxYspP128vbykuG5DBvVgkXzxxOQZ4GLYqItBfrMwxHzWwrcLtz7okY9y29YMrgCRRnD+ZAzWFWHlhHaR8ay1BxrIFnF+1gwZKdHKsNn01IT/Vz8czhXH7uSMaXDNDcCCIiJxHrhOFyYBJwr5n5nXOPx7h/iTMz49KxF/Dw6id5cdsbfSJh2Lqnkr+9vo03Vu1rG5swemgu750ziotnDiczXZccREQ606OEwczO6OgBU865F4EXzexhoCuTO4kHzR01m0fX/j2pBz8651iz5QiPv7SZNVuPAOGxCeedVcy1F45l8uiBOpsgItINPT3DsNLMHgT+0zl3NHplpOxd5ZIcctNz2g1+XMzHpl6T6JC6LBRyvLX+IH9+aTNb9lQCkJEW4IrZI7nmgjEMGZiZ2ABFRJJUTxOGc4D/BraY2V3AL51z7342rySt44MfF/Hhye/1/OBH5xxLNxzi4ec2smP/MQBys1K59qIxXH3+GN3pICJymnr6LIm1wDwz+wDwI+CLZvY159yCWAYniTNl8ASKsgs5WFPm6cGPzjlWbi7jj89tZPPuSgAG5qbzkUvGc9m5I0hP1Z3DIiKxcFq/TZ1zfzWzZ4GvAo+a2SLg/3Y0vkGSi5kxb+yFPLz6SRZsWejJhGH7vip++/Q6Vm8Jj1EYkJ3GRy8dz5VzRpGa4u0zIiIiySYWf35lAsuB3wM3AWvM7FfAdzQnQ3K7ZMx5/HndM6w99A57qw4wPK840SEBUF5Vzx8WbOTlZXtwDrIzUvjwJeO55vzRpKfpjIKISDz09C6JW4BZkddYoAlYBdwfef80sMHMPuSceysWgUrvy07N4qJR5/LCttd5dssr3FD6yYTG09Qc5KmFW/nzy1tobAoS8BtXnz+Gj182gZzM1ITGJiLS1/X0z7GvAUuAB4A3geXOuaZ26x8ys28AvwWmnF6IkkjvnfAeXtj2Oq/tfJNPnvl+stOyEhLHincO86un1nDgSC0Qvj3yc1dPZuig7ITEIyLS3/R00GNX5lj4HfCDnvQv3jE8t5hpRZNYfXAjL21fxPsnXd6r2y87Ws///G0tS9YeAMLTN3/xQ2dx5rhBvRqHiEh/F88LvmXAJXHsX3rJe8dfwuqDG3lu60KumXhpr9xiGQo5nn9rF797ej31jS2kp/r5l8vP4NqLxhDw++K+fREROVHcEgbnnANe7W47M7sW+A5QD/iBm51zy05S9zrgm8DBqFVzgE84556K1KskPLaivQecc491N77+aHrxZIpzBnOg+jBL961mdsnMuG7vwJFafv74KtZuC9/9MHtqETd+8CwGDdADoUREEsVTQ8rN7GzgEeAc59wGM7sGeN7MpjjnopOCVvc45+a362MM4bs2nmtXZ5Vzbm6cwu7zfObjvePfw29XPMaCLa/ELWEIhRzPLNrO7/+xkabmIHnZqfz7h87i/LOGahpnEZEE81TCANwOPO+c2wDgnHvGzA4Rvl3zOx3U/2sHZZ8DHnXO1cctyn7o4lGz+dPav7GxbCs7ju5hdH5sHxVScayB+x9dyYpNh8PbmzGcL3xgKnnZaTHdjoiI9IzXLgbPA6IvPywFLuuosnOu0jlX2bps4T9DP0v47gyJoYyUdN4z+jwAnt70Ykz7fmvdAb7841dYsekwOZmpfOu6Wdz66bOVLIiIeIhnEgYzGwjkAQeiVh0ExnSxm/cANc65pVHlRWb2mJm9ZmYvmtmXzExTAXbT1RMuwWc+Fu9exuHa8tPur6k5yANPrObu373Nsdompk8o5Oe3zmXOmUNjEK2IiMSSZxIGoPUG/8ao8kbCs0l2xXXAbzoo3wrc4Zy7CLge+Arws446MLMbzGyZmS0rKyvr4mb7h8KsAi4YMYuQC/H0Oy+cVl+HKur4xv97nWcX7yTg9/Fv107lri/MoSBPAxtFRLzISwlDbeQ9+jx0GlDXWWMzywHeBzwcvc45d41zbkvk807gXsIPzBrYQd0HnXOlzrnSwsLC7u1BP3DtGeGrQy/vWExVw7Ee9bFs4yFu+elCtu6tYvDATH705Qv5wMVj8fk0sFFExKs8kzA45yqASqAoalURsK0LXXwMeMk5d6QLdbcBRtcvdUjEiAHDOHvomTQHm1mw5ZVutQ2FHA8/t5G7/vdNauqbmTV5CPd/9WLGlQyIT7AiIhIznkkYIl4ESqPKSiPlnbmODi5HmNmlkbkd2hsWed/d3QAFPjjpSgCe3/Iqdc1duxmlobGFex5aymMvbMZn8Jn3TuKOz59Ltp4BISKSFLyWMNwDXGFmkwDM7CqgGPhFZPluM1tnZuntG5nZOGAU8HwHfZYAt5lZZqRuHuHHcT/hnDscrx3pyyYMGsOkwvHUNtfz4rbXO61fdrSeb/y/N1iy9gBZ6QH+4wtz+Ni8CboEISKSRDw1D4NzbrmZfYrww6taZ3q8ot2kTemEB0BGf9NcBzzknAt10O1LwEzgFTNrAHKAl4G74rAL/cYHJl3OxrItPLPpJa4c/x5S/Skd1tu8+yh3//YtjlY3Ujwoi+/867mUDMnp5WhFROR0WXgGZ+lIaWmpW7asw1mp+z3nHF9//vvsqtrHDaWfZN7YC99V5611B/jhH5bR1BLirHGD+ObnZukx1CIi3nbSU79euyQhScLM+MDkKwB4auPztARbTlj//Ju7+MH8t2lqCXH5uSO564Y5ShZERJKYEgbpsTnDz2ZYbhFlteW8smMJED7z8NgLm/h/f15FyMG/XD6RL310mp4wKSKS5PRbXHrM5/Px0SnXAPDkhgU0Njfx4FNrefi5dzCD//Phs/jkFWfowVEiIn2ApwY9SvKZXTKDERuGsbtqH3c+8RgbluYS8Pu49dNnc/5ZmuJZRKSv0BkGOS0+8/GRyVcBsL15OWlp8J9fmK1kQUSkj1HCIKclGAzx2qtBQrW5WGojV15jTBuvKbVFRPoaJQzSYy3BED/643JeX7Uf3+GJALx1+A0aWqKfHyYiIslOCYP0SDDk+O8/rWDR6v1kpge46xPvZ+zAkVQ1VvPPra8mOjwREYkxJQzSbc45fv3kGl5buY+MND/fu/E8Jo0q4ONT3weE52WoaartpBcREUkmShik237/jw0sWLKTlICP7/zrbCaMyAdgWtFkpgyeQG1THU9teC7BUYqISCwpYZBu+fNLm3nila34fcY3PzeLM8cNaltnZnxm2ocBWLBlIYdruvKkcRERSQZKGKTLnn9zJw89uxEz+Oq/zOScyUXvqjNm4AguHHkOLaEW/rT2bwmIUkRE4kEJg3TJincO88sn1gDwxQ+dxcUzh5+07ifOvJYUX4BFu5extXxnL0UoIiLxpIRBOrVjfxX3PLSUUMjx0UvH897zRp+yfmFWAVdNuASAP6x+Ej0RVUQk+SlhkFMqr6rnrv99k/rGFi6aPoxPXzmpS+0+OOlKclKz2Fi2heX718Q5ShERiTclDHJSdQ3N3PW/b1Je1cDk0QO5+RMz8Pm69iCpzNQMPjwlPGX0H1Y9SXOwOZ6hiohInClhkA6FQo6fPrKCHfuPMXRQFt/+/Lmkpvi71cflYy9iWE4RB2oO84/NL8cpUhER6Q1KGKRDj724mbfWHyQrI4X/+MJscrNSu91HwB/g8zM/BsAT65/lSF1FrMMUEZFeooRB3uXtDQd55Pl3MIPbPn02Qwdl97ivs4omMXv4TBqDTTy06okYRikiIr1JCYOcYO/han7yx+UAfOa9kzj7jCGn3ednZ3yYNH8qb+5ZwZqDG0+7PxER6X1KGKRNXUMz3//d29Q1tHD+WUP5yCXjY9LvoMyBbQMgf7viMVqCLTHpV0REeo8SBgHCD5T62WOr2Hu4hhFFOdz8iRmYde2OiK64ZsKlFOcMZn/1IQ2AFBFJQkoYBIDnluxk0Zr9ZKQF+PZ155CRFohp/wF/gH+d+XEA/rL+H3rOhIhIklHCIOzYX8X//G0dAF/+6HSGFvZ8kOOpTCuazHklZ9MYbOLBZY9oBkgRkSSihKGfa2hs4Yd/WEZzS4jLzx3JhTOGxXV7n5/5MbJTs1hzaCMLdyyJ67ZERCR2PJcwmNm1ZrbUzF4zs0VmVtpJ/XfMbGHU6+aoOrlmNj/S7wozu9fMYnvOPUn96qk1beMWvvCBqXHfXl56LtfN+CgAD636C0frq+K+TREROX2eShjM7GzgEeBzzrmLgP8Cnjezdz9H+biDzrm5Ua/7o+rMB/zOuVnAbOBC4Ltx2IWksnD5Hl5auofUFD9f/0wp6am9k0NdOPIcZhRPpba5nt+seLRXtikiIqfHUwkDcDvwvHNuA4Bz7hngEHBTTzs0s6nAB4EfRvpsAu4DbjGz+FysTwJlR+t54MnwQ6Fu+MBURhbl9tq2zYwvlP4L6YE03t67ijf3rOi1bYuISM94LWGYByyLKlsKXHaafTYA66L6zAAuOI1+k1b4FsqV1DW0cO6UIi4/d2SvxzAocyCfnvZBAH6z/FGONVT3egwiItJ1nkkYzGwgkAcciFp1EBhziqZZZvbbyJiHhWZ2h5mlt1s/BjjkThySf7Ddun5nwZKdrNpSRk5mKjd9dFpM51vojnljL2Ry4XiqGqv59bI/6q4JEREP80zCAGRF3hujyhuBzFO02wQ8EBnz8BHgKuAvUf121Ccd9WtmN5jZMjNbVlZW1tXYk8aBI7X89un1ANz0kWnk56R30iJ+fObjpnM/R0ZKOkv3reaVHYsTFouIiJyalxKG2sh7WlR5GlB3skbOuU8755ZGPh8BvgNcbWbT2/XbUZ901K9z7kHnXKlzrrSwsLB7e+BxwZDjvkdX0NgU5KIZwzh/2tBEh0RhVgHXz/wXAH638s8crD6c4IhERKQjnkkYnHMVQCUQfUdEEbCtG1211h0Xed8ODLYTz7sXRdXtF/7+2jY27KggPyeNf//QWYkOp80FI2dx3ohSGlsa+fmbvyMYCiY6JBERieKZhCHiRSB63oXSSPm7mNmZZnZ9VHHrzEO7I+8vEB7gOCWqz3pg0WlFm0QOltfy8ILwkyK//LHp5GSmJjii48yM68/+BAUZ+Wyp2MmTGxYkOiQREYnitYThHuAKM5sEYGZXAcXALyLLd5vZunaDGguAr5tZQWR9GuFbM98ClgM459YDTwG3ReqkADcD9znnanprxxLJOccv/7KappYQc88ezqzJp5rWIjGyU7O46dzPYRhPbFjAhsNbEh2SiIi046mEwTm3HPgU8JCZvQZ8G7jCOdd6V0M64YGKrZcX1hAe4PismS0kfMbgIHCtc679ee3rAMxsKeFkYjFwZ1x3xkNeXbmPlZvLyMlM4d/eF//ZHHtq6pCJvH/S5YRciPuX/IbKhmOJDklERCJMt7KdXGlpqVu2LHpaiORSXdfE/7n3ZSprGvnKx6ZzWQLmXOiOYCjIdxfez8ayLUwdPJE7Lv4KPp+n8loRkb7spPfZ6zdxHzf/mQ1U1jQyZUwB884ZkehwOuX3+bllzr+Rl5bDusObeHz9M4kOSUREUMLQp63fXs4/39pFwG/c9JHETdDUXfkZedw8598wM57csIBVB9YnOiQRkX5PCUMf1RIM8Yu/rALgI5dMoGRITmID6qapQyby8anvA+Bnb/6Ow7XlCY5IRKR/U8LQRz27aAd7DtVQPCiLj146PtHh9MgHJl3BjOKp1DTV8sPXH6ChuSHRIYmI9FtKGPqgqppGHnn+HQCuv3YqqSn+BEfUMz7z8ZXZn2dozhB2V+3jZ2/NJ+RCiQ5LRKRfUsLQB/3xuXeobWhhxoRCZk0ekuhwTktWaiZfv/CLZKVksGzfah5b+3SiQxIR6ZeUMPQxO/ZX8fybO/H5jOvfPzVpBjqeytCcIXz1vC/gMx9PbXyON3a9neiQRET6HSUMfYhzjv/56zpCDq4+fzQjinITHVLMnFU0ietmfBSAB97+A5uPbE9wRCIi/YsShj5k8doDrN12hJzMVD55+cREhxNzV4y7mHljL6Q51MK9r/+S/ccOdt5IRERiQglDH9HUHOS3T4fnK/j0e88g20MPl4oVM+PfZn6cmUPPpLqplu+/+nMq6isTHZaISL+ghKGPeHbxTg5X1DGyKIcrPD798+lonQly/MBRlNVV8F+v/YK6pvpEhyUi0ucpYegD6hqaefzFzQB89urJ+P19+7CmB9L4xkU3UZwzmF2Ve/nxol/THGxOdFgiIn1a3/5m6SeeWriN6romJo0ayKxJyX0bZVflpmXz7Yu/woD0XNYd3sR9S35DSyjYeUMREekRJQxJ7mh1A399dSsAn7t6cp+4jbKrBmcV8K2LvkxWaiZL963m52/+jqCSBhGRuFDCkOQef3EzDU1BZk0ewpQxBYkOp9eNyh/Oty/6Mhkp6SzZs5xfvv0QoZBmgxQRiTUlDEnsYHktzy3ZiRl89qrJiQ4nYcYVjOJbF32JtEAar+96m18v+6OmkBYRiTElDEnsj8+/Q0vQ8Z6zSxhV3HcmaeqJiYPGcvuFN5HqT+GVHYt5cOkfdaZBRCSGlDAkqV0Hj/Hqir0E/D4+ecUZiQ7HEyYPHs83Lvw/pPpTeHnHYn721u80EFJEJEaUMCSpx1/YjHNw+bkjGDIwM9HheMaZQ87gWxd9mYxAOot3L+Mni35Nk265FBE5bUoYktCeQ9W8vnofAb/xkUsmJDocz5k8eDzfmXsz2alZLN+/lntf/wUNzQ2JDktEJKkpYUhCj78UPrsw75yRFOZnJDocTxpXMIr/fM9XyUvPZe2hTdy18D4qG44lOiwRkaSlhCHJ7C+r4bUVe/H7jI9cMj7R4XjaiAHD+O4lX2NwVgHbKnZxx4s/1AOrRER6SAlDknn8pc2EHFxSWqKxC11QnDOYu+d9nbH5IzlcW84dL/2Yd8q2JjosEZGko4QhiRwsr+WV5Xvx+YyPXqqxC101ID2X/7jkq8wceiY1TbV8b+H9LNq9NNFhiYgkFSUMSeTPL20hFHLMnTmc4kFZiQ4nqaQH0rjt/Bu5fOxFNIdauH/Jb/nj6qc0V4OISBd5LmEws2vNbKmZvWZmi8ys9BR1h5jZD8zsDTNbaGYrzex2MwtE1auMrG//+nj89yZ2Dh+t46Wlu/EZfGyezi70hN/n59/O/gTXzfgoPvPxt3f+yb1v/JLaprpEhyYi4nmBzqv0HjM7G3gEOMc5t8HMrgGeN7MpzrmORqt9FLgCmOucqzazEmAFkAHc2a7eKufc3DiHH1d/f207wZDjohnDGFaYnehwkpaZcdWESxiRN4z/Xvw/rDywnm+9cC+3XfDvDM8rTnR4IiKe5bUzDLcDzzvnNgA4554BDgE3naT+YeBHzrnqSP09wOPAp3oh1l5TU9/MP9/aCcCH5o5LbDB9xNQhE/mvy29n5IDhHKg5zO0v3MOrO95MdFgiIp7ltYRhHrAsqmwpcFlHlZ1zjzvnHo0qrgfS4hBbwjy3ZCf1jUGmjR/E2OEDEh1OnzE4q4C7L72NC0eeQ2OwiV+8/Xt++dZDNLQ0Jjo0ERHP8UzCYGYDgTzgQNSqg8CYbnQ1h/BZhvaKzOyxyLiIF83sS2bmP41we01zS4inX98GwAd1diHm0gKpfOnc6/j3WZ8h1Z/Cwp1LuP2Fe9hduS/RoYmIeIpnEgagddh/9J93jUCXJhwws3nACOC7Uau2Anc45y4Crge+AvzsJH3cYGbLzGxZWVlZV2OPm9dW7qXiWCMji3KYOXFwosPpk8yMS8acxw/mfYNhuUXsO3aQ21+4h2c2vaTHZIuIRHgpYaiNvEdfTkgDOh3GbmYjgQeAa51zle3XOeeucc5tiXzeCdwLfDFyVoOoug8650qdc6WFhYXd3olYcs7x1MLwJEMfnDsOM0toPH3diAHD+K/Lvsklo8+jOdTCQ6v+wndfuY/DteWJDk1EJOE8kzA45yqASqAoalURsO1Ubc1sCPA34Hrn3MoubG4bYHTvUkevW7HpMLsOVjMwN52LZgxPdDj9QnogjX8/5zN8/YIvkpeWw4ayLdz23N28vH0RzrlEhycikjCeSRgiXgSi510ojZR3yMzygWeA251zr0bKbmi3/lIzuzaq2bDI++7TjjiOWs8uvO/CMaQEvHao+rbSYWfxkyu/wznDp1Pf0sCvlj7M9xbez4Hqw4kOTUQkIbz2LXQPcIWZTQIws6uAYuAXkeW7zWydmaVHlrOBBYQTijIzK41M9HRjuz5LgNvMLDPSJg/4KvCEc86zv/2376ti9ZYjZKT5uXLOqESH0y/lpufwtfNu4Mvnfp6ctGzWHd7Erc99jyc3LKAlFEx0eCIivcpTEzc555ab2aeAh8ysHvADV7SbtCmd8ADI1ov5NwPnRl7fPEm3LwEzgVfMrAHIAV4G7orPXsTGM29sB8KPsM7OSElwNP2XmXHhqHOYVjyZP6x6gld3vsmja//Oot3L+PyMjzF1yMREhygi0itM12VPrrS01C1bFj0tRPxV1zVx3V3P09QS4tffvJShmtnRM9YeeocHlz3CoZrwHTSzS2bymWkfojCrIMGRiYjExElH13vtkoQAL769m6aWEDMmFCpZ8Jgzh5zBT664g49PfR9p/lTe3LOCWxbcxePrnqGhuSHR4YmIxI0SBo8JhRzPLt4BwNXnj05wNNKR1EAqH55yFf991X9w3ohSmoPN/GX9P/jys//BP7e+qvENItInKWHwmBWbDnOwvI7B+RmUTo6+w1S8ZFDmQG6Z82/cdcn/ZdzAUVQ1HON/lz/K1xZ8lyV7lmvSJxHpUzw16FHgH4vCZxfee95o/D5N1JQMJhWO5/vzvs5be1fypzV/40DNYf578f8yMm8YH55yFecMn47PlJuLSHJTwuAhB8trWf7OIVICPi47Z0Siw5FuMDNml8ykdNg0Xt6+iCc3LGBX1T5+uvh/GJE3jA9PeS/nDp+hxEFEkpYSBg95dvFOnIMLpw8jL7tPPXCz3wj4/Fw+7iLmjp7DK9sX89eNz7O7ah//vfh/Kc4ZzPsmzuOiUbNJ9etWWRFJLrqt8hR687bKxuYgn//u81TXNfOTmy9iwoj8XtmuxFdzsJlXdizmbxv/SVldBQB5aTlcOX4ul427iNw03QUjIp5y0mvhOsPgEa+v3Ed1XTPjSwYoWehDUvwpXD7uYi4dcwFv7l3B3ze+wI7KPTy27mme3LCA80fO4spxcxkzUJegRMTblDB4xD/f2gXAezUNdJ/k9/k5f8QszispZe2hd/jH5pdYeWA9C3csYeGOJUwoGMNlYy9kdslM0gKpiQ5XRORdlDB4wJ5D1WzcWUFGmp8Lpg/rvIEkLTPjrKJJnFU0iYPVh3l+62u8smMxm8u3s7l8O79d+RgXjjiHS8acr7MOIuIpGsNwCr01huF3T6/nyYVbueycEXzl4zPivj3xloaWRhbtWsrL2xexpWJnW/nIvGFcOOoczh8xi4JMXaYSkV5x0jEMShhOoTcShpZgiM9/959U1jTyoy9fyBmjBsZ1e+Jtuyv38fKOxby28y1qmmoBMIwpgydw/ohSzhk+nRwNlBSR+FHC0BO9kTAsWXuAH8x/m5Ih2fzitksw02RNAi3BFlYeXM/rO99m+f41NIdaAPCZj6mDJzK7ZCbnDJtGbnpOgiMVkT5Gd0l41Ytv7wbgsnNGKlmQNgF/gFnDpjFr2DRqm+p4a+9KluxZztpDm1hzaCNrDm3kf5Y/wsSCMZQOO4vSYdMYmjMk0WGLSB+mMwynEO8zDBXHGvj89/6JAfPvvIIBOZqsSU6turGGpftWs2TPCtYd3kSw3YOuirMHM714CtOLJzOlcAKputtCRLpPZxi86KWluwmFHHPOLFayIF2Sk5bNJWPO55Ix51PXXM/qgxtYum8NK/ev5UDNYQ5sOcyCLa+Q4k9hcuE4pg4+gzOHTGTUgBJ8Pk1LLSI9p4QhQZxz7S5H6PY56b7MlAzmlJzNnJKzCYaCbCnfwaqD61l1YAPbj+5m9cGNrD64EYCs1EwmF45nUuF4JheOY+SA4fh9/gTvgYgkEyUMCbJ+ezn7j9QyMDedmRMHJzocSXJ+n58zCsdxRuE4PnHm+6lqOMbaQ5tYd+gd1h7eRFltOUv3rWbpvtUAZATSmTBoNBMKxjBh0BjGDxxNZmpGgvdCRLxMCUOCvLR0DwCXzirB79epYomtvPRcLhg5iwtGzgLgUE0ZGw5vYWPZVjYe2cqhmrITzkAYxrDcIsYNHMWYgSMYO3AkIwcM10OyRKSNEoYEaGoOsnjtfgAuKS1JcDTSHwzJLmRIdiHvGXMeABX1lWw6so3NR3awuXw7O47uYe+xA+w9doCFO5cA4Dcfw3OLGZk/nFEDShidX8KIvKGaB0Kkn1LCkABLNxyirqGFccPzGD5Y99FL7xuYMaBt/ANAU7CZnUf3sP3obrZV7GJ7xS72Vh9kV9U+dlXt4zXeamubn55HSd5QRuQNZXheMcNyixiWW0R2alaidkdEeoEShgRYuCJ8OeLimTq7IN6Q6k9hwqDweIZWDS2N7K7cx87Kvew8uoedlXvZU7Wfow1VHG2oYs2hjSf0kZeey9CcIRTnDGZozmCKc4ZQlF3IkKxBusVTpA9QwtDLquuaWLbxMD6Di2boQVPiXemBtHclESEXoqy2nN1V+9lTtZ+9xw6yr+oA+6oPUtVwjKqGY2ws2/KuvvIz8ijKLmRw1iAGZxVQmFXQ9p6fMYCA7tgQ8TwlDL1s0er9tARDTB9fyMDc9ESHI9ItPvO1jYeYNWxaW3nIhThSd5QD1Yc4UH048jrEwZoyymrLOVpfxdH6KjaWbX1Xn2ZGfnoegzIHUpCZT0HGAAZmDmBgRj4DMwaQn5FLfnqezlKIJJgShl62cMVeAC6eOTzBkYjEjs98DI6cNZhWNPmEdcFQkPK6oxysKeNwbTll7V6H68qprD9GRX0lFfWVUH7ybWSlZDAgI48B6bnkpeWQl54b/pyeQ25aDrlp2eSl55CTlk1GIF1TrYvEmOcSBjO7FvgOUA/4gZudcyedn9nMcoGfAVMi9V8Avu2ca2lXpxj4FTAYSAMecc79OG47cRKHj9axfns5qQEf551V3NubF0kIv8/P4OxBDM4e1OH6llCQivpKyusqOFJ7tC15KK8/ytG6So42HONoQxW1zfXUNtez79jBLm0zNzWb7LQsclKzyE7LIju19ZVJVkomWamZZKdmkpmSQVZqJlkpGWSmZBDwe+7XoogneOonw8zOBh4BznHObTCza4DnzWyKc+5kvyXmA7XOuVlmlgosBL4LfCvSpw94GljgnPuOmeUBK8zsmHPuwTjv0gleW7kPgHOmFJGZrvvbRQACPn/b2QkKO64TciFqmuqorK+isuEYlQ3HqGqoprKhimONNeFXQzXHGqs51lRLY0tj2+DM7krxp5CZkkFmIJ2MlMgrcPw9PSWd9EAaGYE00gNppKekkeZPIz2QSlqkLM2fSmoglXR/Kqn+VE3LLX2CpxIG4HbgeefcBgDn3DNmdgi4ifBZhxOY2VTgg8BZkfpNZnYfMN/MfuCcqwGuAqYD8yJ1qszs18AdZvY/rhefvvVq5HLEXF2OEOkWn/nITcsmNy2bEXQ+WLgp2ExNYy3HGmuoaapte1U31lLbXEdNUx21ra/mOuqaG6hrqqO2uZ7mYDNVwWaqOBaz+FN8AVL9KaQGUsPJhD+VFH+AVH9quNyfQoo/hVRfCin+ACm+ACn+1s/h94DveHnAF4iU+dvK/ZHPrWV+n5+A+Qn4/Pgjr4CF33W5RnrCawnDPODeqLKlwGV0kDBE6jcA66LqZwAXAM9F6mxzzlVG1SkBJgLvxCLwzuw8cIydB46Rk5nCzDP0GGKReEr1p4QHTmYO6FY75xxNwWbqmuupb64PJxLN9TS0NFLf3EB9SwMNLY00tDTQ0NwYLm9ppLGlkcZgEw0tjTS2NLUtN7Y00RhsojnUQnOohdrm+vjscDf5zNeWQPh8vqh3P37z4Tffuz77Wj+bD7/Ph7V+jrz7zCLv7V8nllnbsmFYh2UnWz7+ubU8XKd1n8Jl4YG0hi/8OdLm+DuARdpF6rerC5yiXvRn2vqObtfx8vF20W1pjfX4yhOXWzuI9Ha8nTEwPa9XLqV5JmEws4FAHnAgatVB4L0naTYGOBR1luBgu3Wt7x312bquVxKGhcvDcy+cP20YKQGdnhTxIjMjLZBKWiCV/Iy8mPTpnKM52ExTsJnGYBNNwWaa2t7Dn5uDLTQFm9vqNYdaaI56bwm2tCUeLZHlltblYAstLkhLKBheFwoSjLxaQuF1rctBFyLkQoSCIZppjsk+SmL96IpvM3JA/M9ceyZhAFqniWuMKm8EMk/RpqP6tGvTlTptzOwG4AaAESNi9xTJyppGfKbLESL9jZmRGgiPacgm8bNhOufaEofw+/FEIuhChEJBWlyQUCicWLTWC7n2y5GkwwUJOUfQRdaHXFsd50KEnGtrF3IhHK3L4XfnXLuyEM5xfDkUwkG4n0gZDkI4XGv7yGfnHCEctH+P9N2+Tvg/wLUrD/+rEHIOcG0xtP5btcZwvF1k/QmfadtWpOkJ2wovH+/3hLaRle2ia7d8woE7sd92Zb315FkvJQy1kfe0qPI0oO4UbTqqT7s2tUBuJ3XaRAZCPghQWloas/ENt3xiJp+9ajIDsqPDFRHpPWZGwB/w1C9/SQ6eOTfunKsAKoGiqFVFwLaTNNsODLYTR/C0tt/Wrk5Hfbav0ysG5qbj82mwkYiIJB/PJAwRLwKlUWWlkfKOvEB4gOOUqPr1wKJ2dcaZ2YCoOnucc5tON2AREZH+wGsJwz3AFWY2CcDMrgKKgV9Elu82s3Vmlg7gnFsPPAXcFlmfAtwM3Be5pRJgAbAK+GqkTi7hMQp399I+iYiIJD1PXcZyzi03s08BD5lZ60yPV7SbtCmd8EDF9uf1rwN+bmZLI/VfBO5s12coMnvkr8xsSaSPB3t70iYREZFkZr04b1HSKS0tdcuWnXRWahERkb7mpAPtvHZJQkRERDxICYOIiIh0SgmDiIiIdEoJg4iIiHRKCYOIiIh0SgmDiIiIdEoJg4iIiHRKCYOIiIh0ShM3nYKZlQG7YtjlIOBIDPtLJO2LN/WVfekr+wHaFy/qK/sBsd+XI865KztaoYShF5nZMudc9MO1kpL2xZv6yr70lf0A7YsX9ZX9gN7dF12SEBERkU4pYRAREZFOKWHoXX3pCZnaF2/qK/vSV/YDtC9e1Ff2A3pxXzSGQURERDqlMwwiIiLSqUCiA+hLzOxa4DtAPeAHbnbOLTtF/VzgZ8CUSP0XgG8751p6IdxT6sG+vAMcjCp+yjl3f/yi7JyZpQLfBW4FxjnndnZS38vHpLv74rljYmbXAP8HSAOygKPAN5xzazpp57njchr74sXjchFwCzCQ8L/vAOB/TxWTmaUA9wIXAw5YCdzinKuNd7yniKkn+7Gwg+Ilzrnb4xBij5jZl4CfA+9xzi08Rb34/pw45/SKwQs4G6gBJkeWrwHKgaJTtHkS+EPkcyqwGPhBku7LwkTH3UFMo4AlwO8J/0Ib1YU2Xj0mPdkXLx6TI8An2y3fA5QBQ5LtuJzGvnjxuPwKuLPd8jQgCFxzijY/BV4h/IenAX8GHknC/fDc8YiKbyjh+YAcMLeTunH9OUn4P0ZfeQF/AZ6IKtsAfO8k9adG/gc4s13Zx4A6IDuZ9iWyfmGij8FJ/o3HAXO78iXr8WPSrX3x8DF5Mmq5MLI/n0m249KTffHwcZkM5ESVlRM+Y9BR/XygCXhfu7JzIvs/Lln2w6vHIyq+J4AbO0sYeuPnRGMYYmceEH3Kfilw2SnqNwDroupnABfEPLru6e6+eJJzbp1zbms3mnj2mPRgXzzJOfehqKL6yHvaKZp58rj0cF88yTm3wTlXDWBmPjP7AtBI+KxBRy4GUjjx98RKwn/Nz4tnrKfSg/3wNDN7H9AMPN+F6nH/OdEYhhgws4FAHnAgatVB4L0naTYGOOQiaWC7+q3rEqKH+wKQZWa/JfxXcAh4Efixc64hLoHGhyePyWlIhmMyh/Avub+fok6yHJeu7At4+LiY2R3Al4AK4Crn3L6TVB1D+K/ZtrEYzrlmMyvHA8ekG/vRWv9+YAbhSytvAnc756riHuipY8oCvg9cQdeS0Lj/nOgMQ2xkRd4bo8obgcxTtOmoPqdo0xt6si8Am4AHnHMXAR8BriJ8aSOZePWY9JSnj4mZGeGBtXc45w6foqrnj0s39gU8fFycc3cDxYS/qF41s/NOUjULaI76coLOf0/0im7sB8BqYEHkeFxJOJFbGBnUmUjfA37lnIv+4+1k4v5zooQhNlpHBUdngWmErx+drE1H9TlFm97Qk33BOfdp59zSyOcjhH95Xm1m0+MRZJx49Zj0SBIckx8Au5xzP+mkXjIcl67ui+ePiwv7I/Aa4YGcHakFUiKJUnun/D3Rm7q4HzjnbnbOPRf5XEv4LqTpwNW9EWdHzGwmcC7hQZxdFfefEyUMMeCcqwAqgaKoVUXAtpM02w4MjvqBa21/sjZx18N96Uhr3XExCKu3ePKYxJBnjomZ3UJ4gNrnu1Dd08elm/vSkYQfl8gtu9E2EL49ryPbCZ++H9KujwBQQAKPSQ/2oyM7iAzejElQPXM14bEHL0du+3w0Un6fmS00s4kdtIn7z4kShth5EYh+YlhppLwjLxD+H6L9/8ilhAdOLYp5dN3TrX0xszPN7Pqo4mGR990xji2evHxMusXLxyQS11XAx5xzLWY2xsxONVDOs8elu/vi4eOyvIOzBUOBk137f5XwXRLtf0/MIHzv/8l+5/WGbu2HmQ02s29HFRcTToYSdjycc99zzs10zs11zs0FPhFZdUukbFMHzeL/c5LoW0b6yovw3AXVwKTI8lWEB9wURZbvJjx6Nb1dmyeB30c+pwBv4I17/ru1L4Rv9dsMFESW04BnCA8e8ntgf+bSwa2IyXRMursvXj0mhH/x7SA8yr408roR+M9kOy492RcPH5edwJfbLZ9N+Ivm5sjyjYTnAihuV+enwEscn4fhMRI/D0O39oPw/CblrT9PhP+I/nXkuCb0Vuqo/RpF1G2Vifg50V0SMeKcW25mnwIeMrPW2RGvcM61jlJNJzzwpH32ex3wczNbyvHM/M7ei7pjPdiXNYQHbT0bqZ8NrAL+1TkX7NXg24mcnvwn4dneAB41s/3u+O1wSXNMerAvnjwmwB8If8EsjCq/q93nZDkuPdkXrx6XbwHXm9knCd8amQF8DXggsj41UuZv1+abhGd6fIvwl9kq4OZeivdkursfB4GfAH8yswbCAwd3AZc652p6M/CTMbP7gNmRxfvMbItz7qMk4OdED58SERGRTmkMg4iIiHRKCYOIiIh0SgmDiIiIdEoJg4iIiHRKCYOIiIh0SgmDiIiIdEoJg4iIiHRKCYOIiIh0SgmDiHiOmX3fzH6W6DhE5DglDCLSa8zsR2b2fBeqXg08G+94RKTrlDCISG86B3j7VBXMbCjhRwsv7I2ARKRr9CwJEYm7yMOzagg/Qa/VRufc5A7qXg98wDl3zSn6+yTw1Uif6cBC59ztsY1aRNpTwiAicWdmPmAGsAw4F9gNNDrnjnZQ90ngBefcA9HrIuuHRtpPcM5tN7NC4B3nXEHcdkBEdElCROLPORcCioFqYKlz7uBJkoVUYB6nHr8whPCje0dF+i4Drop1zCJyIiUMItJbZgCr3alPa14I7HHO7TpFnVXAH4AXzexlM7sBWBe7MEWkI0oYRKS3TAdWdlLnKjq5O8KFfRY4E1gOfB9YZWYDYhCjiJyEEgYR6S3TgDWd1LkK+MepKpjZMDOb45xb75y7DZgCDAUujU2YItIRJQwi0lsCwBlmNrSjswFmNprwF/8bnfQzHrjXzAKRZR9gwJYYxioiUZQwiEhv+TbwCWAv8F8drL+a8N0RLZ308w6wDVhiZguBp4GbnHOdnb0QkdOg2ypFxBPM7FngL8653yY6FhF5N51hEBGvWEj4bIGIeJDOMIiIiEindIZBREREOqWEQURERDqlhEFEREQ6pYRBREREOqWEQURERDqlhEFEREQ6pYRBREREOvX/ATkvgi2C/BWHAAAAAElFTkSuQmCC\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 backward Euler scheme. The backward Euler scheme proceeds fully implicitly with the *recursive* update formula\n",
    "\n",
    "$y_{n+1}^{i+1} = y_n + f(y_{n+1}^i) (t_{n+1} - t_n)$\n",
    "\n",
    "We define the backward Euler update as\n",
    "\n",
    "$\\Delta y_\\text{BE}^{i} = y_{n+1}^{i+1} - y_n = f(y_{n+1}^i) \\Delta t$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Backward_Euler_Update(y_i,dt):\n",
    "    return dy_dt(y_i)*dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks exactly like what we had for the forward Euler. The first backward Euler step will in fact be exactly the same as a forward Euler step. But we do not stop there.\n",
    "\n",
    "We now need to define the time stepping itself. So we choose a time-step size and the same starting value (initial condition) of $y_0 = y(t=0) = 0$ to initialize our solution field."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, for every time step we append to this array the value of the next time step by applying the backward Euler update rule. Each new time step value then serves as the initial condition for the following time step.\n",
    "\n",
    "We will *recursively* re-fresh the Backward-Euler update with a new estimate for the rate until the estimated solution does not change anymore, i.e.\n",
    "\n",
    "$\\left| y_{n+1}^{i+1} - y_{n+1}^{i} \\right| < \\epsilon_\\text{abs} \\cup \\left| \\frac{y_{n+1}^{i+1} - y_{n+1}^{i}}{y_{n+1}^{i}} \\right| < \\epsilon_\\text{rel}$\n",
    "\n",
    "with suitably chosen absolute and/or relative tolerances.\n",
    "\n",
    "For simplicity, let's set only an absolute tolerance of $10^{-6}$ here.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrate(dt):\n",
    "    #Startwerte\n",
    "    t_end = 4.\n",
    "    absolute_tolerance = 1.e-6\n",
    "    max_iter = 1000\n",
    "    iteration_counter = np.array([0])\n",
    "    y = np.array([0.])\n",
    "    times = np.array([0.])\n",
    "    #\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_old = y[-1] #Starting value for recursive update\n",
    "        i = 0\n",
    "        #\n",
    "        while True:\n",
    "            y_new = y[-1]+Backward_Euler_Update(y_old,dt) #one backward Euler step\n",
    "            i += 1\n",
    "            if (np.abs(y_new - y_old) < absolute_tolerance or i > max_iter): #if change is below tolerance, stop iterations\n",
    "                break\n",
    "            y_old = y_new #preparation of next recursion\n",
    "\n",
    "        y = np.append(y,y_new) #append the new found solution to the solution vector\n",
    "        iteration_counter = np.append(iteration_counter,i) #store how much iterations this time step took to converge\n",
    "    return times, y,iteration_counter\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "00814c34bdac420c8f4ef209c08e04de",
       "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=1e-3,max=1,step=0.1,description='$\\Delta t$ / s'))\n",
    "\n",
    "def plot_BE(dt=0.1):\n",
    "    time_n, y_n,iter = integrate(dt)\n",
    "    fig, ax = plt.subplots(ncols=2,nrows=2,figsize=(18,12))\n",
    "    #solution\n",
    "    ax[0][0].plot(time,y_analytical(time,0),label='analytical')\n",
    "    ax[0][0].plot(time_n,y_n,label='Backward Euler',marker='d',ls=':')\n",
    "    ax[0][0].legend()\n",
    "    ax[0][0].set_xlabel('$t$ / s')\n",
    "    ax[0][0].set_ylabel('$y$')\n",
    "    #error\n",
    "    ax[0][1].plot(time_n,y_n-y_analytical(time_n,0),ls=':',marker='d',label='absolute')\n",
    "    ax[0][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[0][1].legend()\n",
    "    ax[0][1].set_xlabel('$t$ / s')\n",
    "    ax[0][1].set_ylabel('error')\n",
    "    #iterations\n",
    "    ax[1][0].plot(iter)\n",
    "    ax[1][0].set_xlabel('time step')\n",
    "    ax[1][0].set_ylabel('number of iterations')\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.2$?\n",
    "- What happens for $k=4$ and $\\Delta t = 0.3$? How can this be fixed?\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
}