03_Generalized_Midpoint_Picard.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "<div style=\"display: flex; align-items: center;\">\n",
    "    <img src=\"https://github.com/nagelt/Teaching_Scripts/raw/9d9e29ecca4b04eaf7397938eacbf116d37ddc93/Images/TUBAF_Logo_blau.png\" width=\"500\" height=\"auto\" height=\"auto\" style=\"margin-right: 100px;\" />\n",
    "    <div>\n",
    "        <p><strong>Prof. Dr. Thomas Nagel</strong></p>\n",
    "        <p>Chair of Soil Mechanics and Foundation Engineering<br>Geotechnical Institute<br>Technische Universität Bergakademie Freiberg.</p>\n",
    "        <p><a href=\"https://tu-freiberg.de/en/soilmechanics\">https://tu-freiberg.de/en/soilmechanics</a></p>\n",
    "    </div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 3 - Generalized midpoint rule, 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",
    "%run plot_functions/plot_settings.py"
   ]
  },
  {
   "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": "iVBORw0KGgoAAAANSUhEUgAAAgwAAAF9CAYAAABoEDD2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABOsklEQVR4nO3dd3xc1Z338c9vZtSbZVm2JfeCjW0MNpYBUw2YEiCkV1LILgvJkgIP8OwmlN1kSd9kSfIkZNlNAiQhQCAhoYMB23TcC+69ypasZvUy5/ljRvJ4LKt5RnNH+r5fr3ndueeec+7vco30073nnmvOOURERES64kt0ACIiIuJ9ShhERESkW0oYREREpFtKGERERKRbShhERESkW0oYREREpFuBRAfgZVdeeaV74YUXEh2GiIhIf7ETbdAVhi6Ul5cnOgQRERFPUMIgIiIi3VLCICIiIt1SwiAiIiLdUsIgIiIi3VLCICIiIt3SY5UiIjJoBYNBysvLqaqqoq2tLdHhxF16ejqjR48mJSWl122VMIiIyKC1d+9ezIzx48eTkpKC2QmnIUh6zjkOHz7M3r17mTBhQq/b65aEiIgMWnV1dYwaNYrU1NQBnSwAmBkFBQU0Njb2qb1nEgYzu8bMnjOzV8zsHTN73sxO70G7XDN70MyWmtkKM/uhmQWi6hSZ2d/M7O1wndvjdyQiIpJMfD7P/CqMu5NJirz0X+lB4A/OuUudc+cAq4FXzGxED9r5nXNzgXOAC4DvtG80Mx/wNLDGOTcPuBj4ipndGPtDEBERGZi8lDAscc49ErH+E2AYcPmJGpjZacBHgB8BOOeagfuAW8wsO1ztKmBWuD+cc9XAfwN32UC//iQiIhIjnkkYnHMfjSpqCC/Tumi2AGgE1kWULQUygPMj6mxzzlVF1RkDTO1rvL316vY3+e+lf6Smqba/dikiIhIznkkYOjGPUDLw9y7qTAQOOudcRFlpxLb25YGodtF1OpjZjWa2zMyWlZWV9T7qE3h1+1u8sv0NdlftjVmfIiIi/cWTCUP4VsHdwF3OuUNdVM0CmqLK2tcze1Gng3PuAedciXOupLCwsHeBd2FsXjEAu6r2xaxPERGR/uLJhAH4HrDLOfeTburVcfwti/b1+l7UibuxQ0YBsLt6f3/tUkREBrDRo0fz05/+9JiytWvXkp6ezvr162O+P89N3GRmtwDTgI/3oPp2YLiZWcRtiZHh5baIOldGtYuuE3fj2hMGXWEQEfG0D972t4Ts9+mffKhX9efNm8fSpUuPKbvlllu44YYbmD59eixDAzx2hcHMbiD0VMOnnHOtZjbRzBZ00eRlQgMcZ0SUlRAaMPlmRJ3JZjYkqs4e59ymmAXfjTHhWxK7a/YTDAb7a7ciIjJARScMTz31FCtXruTb3/52XPbnmSsMZvZp4E7gemBm+InHOUARsDBc517gw0CJc67ROfe+mf0VuAP4opmlAN8A7nPOtT+O8DywCrgV+DczywVuBO7tnyMLyU7NoiAzn8P1lZTWlVGc0930EiIikgi9/Us/Uc455xxuu+02KioqyMrK4vbbb+eee+6hoKAgLvvzTMIA/J5QPIuiyiNTpXRCAxUj50+4HviFmS0F/ISSi3vaNzrngmZ2LfBrM3s73McDzrkHYn0A3RmbN4rD9ZXsrtqnhEFERE7KnDlzSE1NZdmyZaxcuZJAIMDNN98ct/15JmFwznX76izn3O3A7VFlNcAXu2m3H7j2pAKMgXFDRrHywDp2V+/jnDFnJjocERFJYmlpacyePZunn36ahx56iEceeaRPb6HsKU+NYRjoxuaFBj7q0UoREYmFefPm8ctf/pJ58+ZxzTXXxHVfShj6UftcDHpSQkREYmHWrFn4fL7jHq+MByUM/ag4dyR+n5+DdeU0tvTt9aIiIiLt/vjHP3LTTTcxY8aM7iufJM+MYRgMAj4/o3OL2FW1l93V+5ky7LiZqUVERLoUDAYpKyvjwQcfZO3atTz22GP9sl9dYehnHbclNOOjiIj0wZIlSygqKuLBBx/kySefJD8/v1/2qysM/WzckFG8vkvjGEREpG/mz5+fkAkAdYWhn7U/KbG7WgmDiIgkDyUM/az9JVS7qvdx7Fu5RUREvEsJQz/LT88jJzWLuuZ6KhqqEh2OiIhIjyhh6GdmFvGqa92WEBGR5KCEIQE046OIiCQbJQwJoEcrRUQk2ShhSIBxQ0YDerRSRESShxKGBBidV4Rh7Ks5QGtba6LDERER6ZYShgRID6QxInsYbS7I/iMHEx2OiIhIt5QwJEj7kxI7q/YmOBIREUlmu3fvJjs7m5UrV8Z1P0oYEmTCkDEA7Kzck+BIREQkmRUXF7Nq1SqmT58e1/3oXRIJMiE/lDDsqFLCICIifRcIBJg8eXLc96MrDAkyPv/oFQZNES0iIl6nhCFB8tPzyEvPpa6lgbK6w4kOR0REpEtKGBLEzJgQno9he+XuBEcjIiLJZvTo0fz0pz89pmzt2rWkp6ezfv36mO9PCUMCTcgfC8BOjWMQEZFemjdvHkuXLj2m7JZbbuGGG26IywBIDXpMoI6Bj3pSQkTEMz752FcSst/HP3V/r+rPmzePX/3qVx3rTz31FCtXruTxxx+PdWiArjAk1HglDCIi0kfnnHMO27Zto6KigqamJm6//XbuueceCgoK4rI/XWFIoBFZw8hMyaCqsYbKhmryM/ISHZKIyKDX27/0E2XOnDmkpqaybNkyVq5cSSAQ4Oabb47b/jyVMJhZKvBt4A5gsnNuZxd1rwf+FSiN2jQP+LRz7q/helXAqqg69zvnHotJ0CfBzBg/ZDTry7awo3KPEgYREemxtLQ0Zs+ezdNPP81DDz3EI488QkpKStz255mEwczGA38CNgP+Hjb7gXPuwYg+JgArgBci6qxyzs2PTZSxNyF/bDhh2M2ZxaclOhwREUki8+bN42c/+xmXXXYZ11xzTVz35ZmEAcgGPg+MBr7Qg/pPdVJ2PfCoc64hdmHFl2Z8FBGRvpo1axY+n++4xyvjwTMJg3NuHYCZje5h/arIdTMzQonGJ2MeXBxNyNc7JUREpG/++Mc/ctNNNzFjxoy478szCUMMXAzUOueWRpWPNLPHgCKgmdCVifudc239HF+ninNGkOJP4VDdYWqb68hOzUp0SCIi4mHBYJCysjIefPBB1q5dy2OP9c+QvIH0WOX1wG86Kd8K3OWcuxC4Afg68PMTdWJmN5rZMjNbVlZWFpdAI/l9fsblhV91XalXXYuISNeWLFlCUVERDz74IE8++ST5+fn9st8BkTCYWQ7wQeAP0ducc9c457aEv+8Efgh8xcyGdtaXc+4B51yJc66ksLAwjlEfpQmcRESkp+bPn08wGGTDhg2ce+65/bbfAZEwEBq38IpzrrwHdbcBBkyMb0g9p4GPIiLidQMlYbieTm5HmNmlZnZtVPGo8NIzb3xqf6fEDr2ESkREPCppEgYzu9fM1plZelT5ZGA88FInzcYAd5hZZrhuHnAr8KRz7lCcQ+6xMXnF+MzH/iMHaWxtSnQ4IiIix/FMwmBmqWa2CLgvXPSomf0loko6kEnodkKk64GHT/DUwyuEZnl8zcwWA68Bi4AvxSruWEj1pzAmtwjnHLur9iU6HBERkeN45rFK51wzML+L7bcDt3dSflcXbfYAX4tFfPE2Pn8Mu6r3sb1yN1OGeWZ4hYjIgOecIzSVz8DnnOtzW89cYRjsJg0dB8D2Co1jEBHpLykpKTQ0JM3kwCetpaWFQKBv1wqUMHhEe8KwtWJnYgMRERlEhg8fzr59+6ivrz+pv76TQTAY5ODBg+Tl9e1Fh565JTHYjRsyGr/Pz76aUhpaGslISe++kYiInJTc3FwA9u/fT0tLS4Kjib+srCyGDRvWp7ZKGDwi1Z/C+LzRbKvcxfbK3cwYPiXRIYmIDAq5ubkdiYOcmG5JeEjHbYnDOxMbiIiISBQlDB4yuWA8ANsqdiU2EBERkShKGDxEAx9FRMSrlDB4yKickaQH0iivr6CqsSbR4YiIiHRQwuAhPp+PieH3Sui2hIiIeIkSBo85Oo5hZ0LjEBERiaSEwWMmDx0P6EkJERHxFiUMHnN04OOuAT/rmIiIJA8lDB4zLHMoeWk51DbXcaiuPNHhiIiIAEoYPMfM9HiliIh4jhIGD+oY+HhYT0qIiIg3KGHwIF1hEBERr1HC4EGTwk9K7KjcQ1uwLbHBiIiIoITBk3LTshmRNYymtmb21hxIdDgiIiJKGLxKb64UEREvUcLgUe0DH7cc3pHYQERERFDC4FlTCiYCsFkJg4iIeIASBo+akD+GFF+AvTUHqG2uS3Q4IiIyyClh8KgUfwoTw+MYdFtCREQSTQmDh00dFrotsal8e4IjERGRwU4Jg4d1jGNQwiAiIgnmqYTBzFLN7Ptm1mpm43tQf6OZLYr6fCOqTpGZ/c3M3jazFWZ2e9wOIMamhK8wbKnYqQmcREQkoQKJDqBdOEH4E7AZ8PewWalzbn4XffqAp4HnnXN3m1kesMLMapxzD5xkyHE3JD2XEdmFHKwtY1fVPiYOHZvokEREZJDy0hWGbODzwO9i2OdVwCzgJwDOuWrgv4G7zMxiuJ+4mdrxeKVuS4iISOJ4JmFwzq1zzm2NcbcLgG3OuaqIsqXAGGBqjPcVF1OHTQJgU/m2BEciIiKDmWcShj7KMrPfmtmS8PiFu8wsPWL7RCD6ZQylEds8r/1JCQ18FBGRREr2hGETcL9z7kLg44RuQTwRsT0LaIpq076e2VmHZnajmS0zs2VlZWWxjrfXRucWkZGSTll9BRX1VYkOR0REBqmkThicc59zzi0Nfy8H7gauNrNZ4Sp1QFpUs/b1+hP0+YBzrsQ5V1JYWBiHqHvH5/MxpWACoHEMIiKSOEmdMHSi/Ub/5PByOzAyqs7IqLqe1z4fw0aNYxARkQRJ2oTBzGaa2Q1RxaPCy93h5cvAZDMbElGnBNjjnNsU5xBjpn3go8YxiIhIoiRNwmBm95rZuohBjQXA/zWzgvD2NOCbwHvA8nCd54FVwK3hOrnAjcC9/Rj6SZtcMB7D2FG5m+bW5kSHIyIig5BnEobwLI+LgPvCRY+a2V8iqqQTGqjYPn/CGkIDHJ8Lt3uT0BMQH3TOtQE454LAtcBsM3sbWAw8kAyTNkXKTMlgbF4xbS7ItspdiQ5HREQGIc/M9Oicawbmd7H9duD2iPUK4Fs96Hc/oaQhqU0ZNpFd1fvYVL6daYWnJDocEREZZDxzhUG61j6OYWNZrOe2EhER6Z4ShiQxrTD04MfG8m0Eg8EERyMiIoONEoYkUZhVQGHmUOpbGthVvS/R4YiIyCCjhCGJTBseGruw/tDmBEciIiKDjRKGJDK9cAoA68u2JDgSEREZbJQwJJHp4SsMG8q2EnQaxyAiIv1HCUMSGZE1jKEZQ6htrmNvdfRLOEVEROJHCUMSMTOmh+dg0G0JERHpT0oYksz0joGPShhERKT/KGFIMu1XGDaUbcE5l+BoRERksFDCkGSKckaQl55LddMR9h0pTXQ4IiIySChhSDLHjGPQbQkREeknShiSkAY+iohIf1PCkIQ65mM4pHEMIiLSP5QwJKHRuUXkpGVT2VhNaW1ZosMREZFBQAlDEjp2HIPeKyEiIvGnhCFJaRyDiIj0JyUMSWrG8NCLqNYd2qRxDCIiEndKGJLUmLxi8tJzqWyoZl+N5mMQEZH4UsKQpMyMmSNOBWDNwQ0JjkZERAY6JQxJ7PRwwrD24MYERyIiIgOdEoYk1n6FYf2hLbQG2xIcjYiIDGRKGJJYQWY+o3JG0tDayNbDOxMdjoiIDGBKGJLczI7bEhrHICIi8aOEIcnNHKlxDCIiEn+eShjMLNXMvm9mrWY2vpu6I8zse2b2hpktMrOVZvZNMwtE1asKb4/8fCquB9KPZhROwWc+thzeQUNLY6LDERGRASrQfZX+EU4Q/gRsBvw9aPIJ4ApgvnPuiJmNAVYAGcA9EfVWOefmxzZa78hMzWDy0PFsPryd9WVbmFM8M9EhiYjIAOSlKwzZwOeB3/Ww/iHgx865IwDOuT3A48B18QnPuzrGMZRqHIOIiMSHZxIG59w659zWXtR/3Dn3aFRxA5AW28i8b6bmYxARkTjzTMIQI/MIXWWINNLMHjOzJWa20My+amY9ueWRNKYUTCAtkMaemgNUNFQlOhwRERmABkzCYGYLgLHAd6I2bQXucs5dCNwAfB34eRf93Ghmy8xsWVlZWdzijaWAP9Dx9sp1BzclOBoRERmIBkTCYGbjgF8B1zrnqiK3Oeeucc5tCX/fCfwQ+IqZDe2sL+fcA865EudcSWFhYXwDjyG9V0JEROIp6RMGMxsO/A34J+fcyh402QYYMDGugfWz9vdKrCndQNAFExyNiIgMNEmdMJhZPvAs8E3n3OJw2Y0R2y81s2ujmo0KL3f3T5T9Y0xeMQUZ+VQ11rCzcm+iwxERkQEmaRIGM7vXzNaZWXp4PRt4HlgIlJlZiZmVADdFNBsD3GFmmeE2ecCtwJPOuUP9ewTxZWbMLpoBwMoD6xIcjYiIDDSeSRjCszwuAu4LFz1qZn+JqJIOZBK6nQDwDeBs4F+BpRGfMyPavAKsAl4zs8XAa8Ai4EvxOIZEm118GgArD7yf4EhERGSg8cxMj865ZmB+F9tvB26PWP8u8N1u+twDfC1GIXrezOFT8fv8bDm8g5qmWnLTshMdkoiIDBCeucIgJy89JZ3phafgcKwpXZ/ocEREZAA5qYTBzM6OVSASG7OLQrclVui2hIiIxNDJXmH4c0yikJg5MzzwcfWB9wkG9XiliIjERrdjGMwseqrljk1Ap5MfSeIU5YxgRNYwDtaVs7ViJ1OGDajpJkREJEF6MuhxAaG3SNZGlRtwYcwjkpMSerzyNF7YuoiVB95XwiAiIjHRk1sSrwO1zrnFUZ9FwJr4hid9MbtY8zGIiEhsdXuFwTn3oS62XRbbcCQWZhROIcWfwvbK3VQ1VDMkIy/RIYmISJLr9aBHMxvVfS1JpNRAKqcNnwrAKj1eKSIiMdCXpySej3kUEnPt00Sv0G0JERGJgb4kDNZ9FUm09oRhdel6WttaExyNiIgku74kDC7mUUjMjcguZExuEQ0tjbxftjnR4YiISJLT1NAD2NzRswBYund1YgMREZGkp4RhADtr1BkALN2/mqDTrI8iItJ3GsMwgE3IH0tBZj6VDdVsq9iV6HBERCSJ9TphcM7NjEcgEntmxtz2qwz7dFtCRET6Lia3JMws18w+YmZKJjym47aExjGIiMhJ6FPCYGbfNbMyM5trZpnAUuD3wNtm9oWYRign5dTCU8hKzWTfkVL21ZQmOhwREUlSfb3CcDEwzTm3FLgOyAfGA5OBm2MTmsRCwOdnTnHowo9uS4iISF/1NWFocM6Vh79/Gvidc67cOVcK1McmNImVs0bNAmDp3lUJjUNERJJXXxOGHDMbZ2YXAxcBDwKYmR/IjFFsEiNnjJxOqj+FLRU7qWioSnQ4IiKShPqaMNwHbAUWAn9wzm0ws3OAVwG9vMBj0gKpnD5yOgDL9umN5CIi0nt9Shicc48AY4E5zrnrw8W7gXuAb8UmNImls/R4pYiInIRAXxs65w4AByLW9wP7YxGUxN6c4pn4zMe6Q5uoba4jOzUr0SGJiEgS0dTQg0ROWjYzhp9CW7BNczKIiEivKWEYROaNKQHgrT3LExyJiIgkG08lDGaWambfN7NWMxvfg/q5ZvagmS01sxVm9kMzC0TVKTKzv5nZ2+E6t8ftADzu7NGz8JuPtQc3UtN4JNHhiIhIEulxwmBm95nZZWaWGo9AwgnCYqAY8Pew2YOA3zk3FzgHuAD4TkSfPuBpYI1zbh6hCae+YmY3xi7y5JGTls3MEacSdEHe1ZwMIiLSC725wlAH/Ag4bGZ/N7ObzGxMDGPJBj4P/K4nlc3sNOAj4ZhwzjUTetzzFjPLDle7CpgF/CRcpxr4b+AuMxuUb908d2z7bYllCY5ERESSSY+fknDO3QncaWZFhH4RXwX8yMx2A88BzwJvOufa+hKIc24dgJmN7mGTBUAjx877sBTIAM4HXgjX2eacq4qqMwaYCmzsS6zJbO6oMwj4Aqwv20JVQzVDMvISHZKIDGDOOZwLL+HY78H2MheuCy78JbLu0W2hCse1idgGx24/pm14Y3R7outHfWlfj6zTWdvj2x27j87+2xyzHt0+qp2j846mjs3H74//CINeP1YZfpzyN8BvwuMFLiSUPPwaKDKzlwklEE875w7HMtgoE4GD7tj/4qUR29qXBzhWZJ1BlzBkpWZyxshpLN+/lnf2ruTKU+YnOiSRpBEMOlrbgrS2BWlpPbpsCzpaw+uhj6MtGFq2tgVpi1hvawvV71gGXcf2yO/B8LZg+NMWsR7aDkF3dN05Qtvd0TbOHa0TdA7X8Z3w9vb6R3+xB8Nlkb/Mg+Ff3B31aK/f3qazxMCd8BelxNZj372KTC8mDJGcc62EZnd8Fbg9PA7hauATwDDgP082wC5kAU1RZe3rmb2oc4zw+IYbAcaOHXvyUXrQuWNKWL5/LW/tXqaEQZJWS2uQpuZWGpraaGxupam5jaaWtvCylaaWYMf35pYgzS1tNLeE6rS0htdbg7S0BGluPVrW0hqkpbV9Gf6EE4NgUL8B+8pngFloiWEGZuFlaNPRckIFx5UbGKEKJ2wTvtncvt7+vaNtuOBo+2P7PyrcNrI/QjEfU6tj+7EVo9tFxtJZ+4hIju3nuPrH303vrzvsJ5UwRHPO7QR+Gf7EWx2QFlXWvl4fUSe3mzrHcM49ADwAUFJSMiB/OpSMOp0Ufwoby7dxuL6Sgsz8RIckg0Qw6GhoaqWuoYW6xhZqG1qob2ihrrGV+sYW6tuXTa00NLbS0BT61De10hj+NDS30djUSluCfnkH/D5SAr7w0vD7faT4fRFLI+APbW//7vdZx3r7d58v9N0f3h79PbQ9sp7hM+tY9/mOXTeLLgdfuNxnR8s61sN1zehYb//e/ovcZ0d/sUd+j95mgPmOTwZ81n+/zCT+Ypow9LPtwHAzs4jbEiPDy20Rda6MahddZ9DJSElndtEM3tu7irf3rOCaqZcmOiRJQo1NrVTXNVNT10R1bTM1daHPkfpmjoSXtfUt1DY0c6Q+nBw0tsTsMrXfZ6Sn+klPC5CW4ic9NUBaqj/0SYlYpvhJ7fj4SEvxk5LiJy3FR0ogtD0l4Ov4pKb4O5KC0Ofo9vZfzCKDUTInDC8D/wXM4OjAxxKgAXgzos5XzWxIxMDHEmCPc25TP8bqOeeOKQklDLuXKWGQDq1tQSpqGqmoaaSyppHKI01U1jRReaSRypomqmubqAp/mpr7NL6ZjLQAWRkpZGekkJWRQmZ6gMy0FDIzAmSGt2WkBchMD5CRdvSTnhYgIzVARnqA9NQAKQFPTSMjMuAlTcJgZvcCHwZKnHONzrn3zeyvwB3AF80sBfgGcJ9zrjbc7HlgFXAr8G9mlktofMK9/R2/15xZfBpp/lS2VOzkYG0ZI7ILEx2SxFlbW5CKmiYOVdZTVllPWVUD5VUNHK5upKyqgYrqRqpqo4f8nFhKwEdeViq52WnkZqUe/WSmkpOVSk5m6JOdmRL6ZKSSlR7ol9HcIhJ7fUoYzOxU51xMnzAITwj1EjAkXPSome13zn00vJ5OaKBi5PXA64FfmNlSQpM9LST0xkwAnHNBM7sW+LWZvR3u44HwOIVBLT2QxtzRs3hj13u8vus9Pj7j6kSHJDFQ29DCgfJaDpTXceBwHQcP13Owop7SinrKqxq6HbTnMxiSk87QvHSG5qSTn5tGfng5JDv0PS8nlSHZaWSkBXR5XmQQsejnQHvUyKyB0MDAf3fOVcY8Ko8oKSlxy5YN3AmOVpeu57uLf8GI7EJ+ftW39cM/SbS2BTlQXsfeQ0fYe6iWfWW17D1Uy/6yOo7UN3fZNj8njeH5mRTmZ1CYn8mwIekMy8tg2JAMCvLSGZKTjt+nfwcig9gJfwD09ZbEWYTGD2wxs28Dv+rrhE2SODOHn0p+Rh4Ha8vYVL6dUwsnJTokiRAMOg5W1LPzQDU799ew6+AR9hw8wv6yWlrbOk/001L9FBVkUTQsi5EFWYwsyGTE0ExGFmRROCSD1JSezrouInKsbhMGMzsDWO+ca2kvc86tBRaY2YeBHxN6P8Ntzrnn4xapxJzP5+OCcWfz940vsWTnO0oYEqilNcju0hq27atm294qtu2rZteBGhpPMLBw+NBMxgzPZvTwHEYNz2Z0YTbFhVkMzU3XlSIRiYueXGF4AziD0COKx3DOPWVmzxEaVPiomb0J/J9Yj2+Q+LlofChheGvPcq6f/QlSA3F5t5hECAYd+8pq2bKnks27q9i8u5Id+2tobQseV3dobjrji3MZPzKXcUW5jB2Rw+jh2aSnJc14ZREZIHryU2c6sL+L7ZnAcuAh4GZgjZn9Grg7/LIn8bAxecVMzB/L9srdLNu/lnPHzkl0SANOU0sbm3dXsmFHBRt2VrBxZwW1DS3H1DGDUYXZTBqdx6RRQ5g0Ko/xxbnkZUfPTSYikhjdJgzOuT3RZWZ2CzA3/JkENBN6fPFn4eV1wHoz+6hz7t3YhSvxcOH4s9leuZslO99RwhADTS1tbNxZwdpt5azbdphNuyqPu3owNDedqePymTI2nyljhzB59BAy01MSFLGISPf6el3zNuBt4H7gHWB5+PXS7R42s38BfktoYiXxsPPHzuX3q55kVel6qhprGJIePZu2dCUYdOzYX82qzWWs2lzG+zsO09J6NEEwg/FFuUyfMJRpEwqYPn4ohfkZGmsgIkmlTwmDc25MD6r9DvheX/qX/pWbnsPsotNYtn8Nb+xaqpkfe6CuoYWVmw+xdP1BVmw8dNyERxOL8zhtcgEzJw1jxsQCcjI1NkREkls8R06VAZfEsX+JoQvHn82y/WtYsvMdJQwnUF7VwDvrDvD22gO8v/3wMS8/GjYkg9lTCpk1pZAzTinU2AMRGXDiljCEXwi1OF79S2zNKZ5JVmomO6v2srNyL+PzRyc6JE8oPVzHG6v389aa/WzZU9VR7vMZMyYWUDJtBHOnjWDsyBzdYhCRAU3PZgkAKf4ULhh7Fi9sXcQr29/gH+d8OtEhJczh6gbeWL2f11fuY9PuoxOZpqb4mXPqcObNLGLutBFk6zaDiAwiShikw4JJ5/PC1kUs2fUu153xEdIDg+eyemNzK++sK+WVpbtZvaWs4xXM6al+zpoxkvPPKGb21OGkp+p/GREZnPTTTzqMHTKKUwomsOXwDt7evZyLJ56b6JDiyjnHlj1VvPjOLl5ftY+GplYAAn4fc6eP4MLZoyiZNkJJgogIMU4YzCxPkzUltwUTz2fL4R0s3P7GgE0Y6htbWLJyH8+/vZPt+47+c50ydgiXzh3LBbNG6akGEZEosf7TqdLMtgLfdM49GeO+pR+cO7aEh1Y9wZbDO9hVtZdxQwbO4Mf95bU888YOFr63u+NqQk5mCpfOHcvlZ49jzIicBEcoIuJdsU4YLgemAT80M79z7vEY9y9xlhZI5YJxZ/Hi1sUs3Jb8gx+dc6zdVs7fl2znvfWlHWMTpk8Yygfmjefc04v1BkcRkR7oU8JgZqd29oIp59xCYKGZ/QHoyeRO4kELJp3Pi1sXs2TXu3zujI+SloQvpAoGHe+sO8CfX9nM1r2h2w4Bv4/5Z47m2gsnMqE4L8ERiogkl75eYVhpZg8A/+6cq4zeGC47rlySw7gho48OftyznPkT5iU6pB5rbQuyeMVennh1C3sP1QKQl53K1edN5Mp548jPSU9whCIiyamvCcNZwH8BW8zs28CvnHNtsQtLEq198OPL215PioShLehYvGIvf3ppI6WH6wEozM/gY/Mns+DscaTptoOIyEnp67sk1gILzOzDwI+Br5jZbc6552MZnCTOvLFzeHDVnz0/+DEYdLy5Zj+PvLix44rCqMIsPrlgChfOHk3A70twhCIiA8NJDXp0zj1lZs8BtwKPmtmbwP/pbHyDJJf0QBoXjjubF7cu5vkti/jy3M8lOqTjrN1Wzm//vq5jjMKIoZl85vKpzD9zNH4lCiIiMRWLpyQygeXAQ8DNwBoz+zVwt+ZkSG4fOGU+L25dzOu73uO60z9MTlp2okMCYO+hIzz4zHrefb8UgKG56Xz68qksmDuWlIASBRGReOjrUxK3AHPDn0lAM7AK+Fl4eR2w3sw+6px7NxaBSv8rzh3JrJHTWVW6nle2v8mHp12R0HjqGlp45KWNPPvGDtqCjvRUPx+/5BQ+dNEkzcYoIhJnff0pexvwNnA/8A6w3DnXHLH9YTP7F+C3wIyTC1ES6aopl7CqdD0vblnMNVMXEPD1/+BB50IDGn/79PtUHmnCZ3DFOeO47opTyc/VUw8iIv2hr4MeezLHwu+A7/Wlf/GO00dOozhnBPuPHOS9vas4d+ycft3/7tIafv2XtazdVg7AtPFD+crHTtc8CiIi/Sye13HLgEvi2L/0A5/5+MApF/ObFY/y/OZX+y1haG0L8uSrW3j05U20tjlys1L50jXTuaRkLD6f9UsMIiJyVNwSBuecAxb3tp2ZXQvcDTQAfuAbzrllJ6h7PfCvQGnUpnnAp51zfw3XqyI0tiLS/c65x3ob32B00fiz+dPav7Hp8Ha2Vexi0tBxcd3f1r1V/PyxlezYXwOEbj988erpeiGUiEgCeWqkmJnNAR4BznLOrTeza4AXzWyGcy46KWj3A+fcgxF9TABWAC9E1FnlnJsfp7AHvPSUdC6ZeB7PbFrI85tf46vnXB+X/bS2BfnTS5t44tUtBIOOEUMz+donZ3HGKYVx2Z+IiPSc155B+ybwonNuPYBz7hngIKHHNTvzVPgT6XrgUedcQ3xCHJyunHwRZsabe5ZR1RD7p2X3l9Vyxy9e5/GFm3HOce0FE/l/t1+sZEFExCO8ljAsAKJvPywFLuussnOuyjlX1b5uZgZ8gdDTGRJDw7OHUVJ8Om3BNl7Yuihm/TrneOndXXz9p4vYuqeK4fkZfP+fz+efPjyT9DRPXQATERnUPJMwmNlQIA84ELWpFJjYw24uBmqdc0ujykea2WNmtsTMFprZV81MLxfopWtPDeVtL25ZTENL40n3V9vQwg8eXsovHl9FU3MbF80ezc9vu5gZEwtOum8REYktL/0JlxVeNkWVNxGaTbInrgd+00n5VuAu59wWMxsPvARMo5NbHWZ2I3AjwNixY3u428Fh6rBJTCuczIayrby87fWOBKIvduyv5vsPLuXA4Toy0wN85aOnM3+O3oguIuJVnrnCANSFl2lR5WlAfXeNzSwH+CDwh+htzrlrnHNbwt93Aj8k9MKsoZ3UfcA5V+KcKyks1P3zaB86NTTb47ObXqGlraVPfSx8bze3/2wJBw7XMbE4j/tuna9kQUTE4zyTMDjnKoAqYGTUppHAth508UngFedceQ/qbgOMnt/qkLDZRTMYlzeKysZqluzs3azfLa1t/OLxVfzssZU0twa57Kyx/OjrF1A0LKv7xiIiklCeSRjCFgIlUWUl4fLuXE8ntyPM7NLw3A6RRoWXu3sb4GBnZnwo/E6Jv298mWAw2KN21bVN3Hn/W7z07i5SAz6+8alZfP1Ts0lL0VASEZFk4LWE4QfAFWY2DcDMrgKKgF+G1+81s3VmdswLBMxsMjCe0NiEaGOAO8wsM1w3j9DruJ90zh2K14EMZPPGnMnwrAIO1B7i3X0ru62/60AN/+dnS9iws4Jheen88GsXsOCs+E7+JCIiseWphME5t5zQmy4fNrMlwJ3AFRGTNqUTGgAZPTfw9cDDzrm2Trp9hdAsj6+Z2WLgNWAR8KVYxz9Y+H3+jgGPT214kdCknp1bur6UO36xhEMV9UwZO4Sf3HIRk0cP6adIRUQkVqyrH/aDXUlJiVu2rNNZqQe95tZmbn7mLqqbjnDnRV/jjJHTj6vz/Fs7+PVf1hB0cOGsUXz907oFISLicSd8WY+nrjBI8kgNpHL11EsBeOL95465yuCc45EXN/KrJ0PJwmcvn8rtn5ujZEFEJIkpYZA+u2LyReSkZrGpfBurSzcA0BZ0/PKJ1fzppU34DL76iVl85opTCU3CKSIiyUoJg/RZRko61556OQCPr3uapuZWfvjwUl58J/QkxDevP4srztHgRhGRgUAJg5yUK065iLy0HLZW7OSbv/8Lb689QFZGCt+56VzOOa0o0eGJiEiMKGGQk5IeSOPqKQsA2G3Lyc1O4fv/fJ7eByEiMsAoYZCT0tDUyjuL0nDNafiyarjuE/lMKM5LdFgiIhJjShikz+obW/j3/3mb97dVkVo5BYDX9r5C0PVs9kcREUkeShikT5pa2viP377L+h2h2Ru/98nPUJCRz67qfby7t/vZH0VEJLkoYZBea20L8sOHl7Ju22GG5qbxvX8+n3Ej8/nI9CsBeHTN32kNdjbppoiIJCslDNIrbUHHfz2ygqXrD5KTGXoaov1tk5dMPI+inOEcqD3Ewm2vJzhSERGJJSUM0mPOOe5/cjVLVu0jIy3Av//TPMaNzO3YHvD5ue70jwDw5/efpb65IVGhiohIjClhkB77wwsbOyZluvsfz2bK2Pzj6swddQanDpvEkaZantr4YgKiFBGReFDCID3y0ru7eHzhZnw+41++OJeZk4Z1Ws/M+PysjwHw7OZXKa+v6M8wRUQkTpQwSLdWbDrEL59YDcBXPno6Z00f2WX9UwomcO6YObS0tfDo2r/3R4giIhJnShikSzsP1PCDh5YSDDo+dvFkrpw3vkftPnP6h/D7/Ly+8z12VO6Jb5AiIhJ3ShjkhA5XN/Dt/3mbhqZWzj+jmC9cNb3HbUdkF3Ll5Pk4HA+veuKY11+LiEjyUcIgnWpuaeO7v3uP8upGpo0fyq2fOROfr3evqP7Y9A+Qk5rF+4c289aeZXGKVERE+oMSBjmOc45fPbmaLXuqGD40kzu/dBapKf5e95OdlsVnT/8wAA+vepKGlsYYRyoiIv1FCYMc57m3dvLK0j2kpvi560tnkZed1ue+Lp54LpOHjqeyoZon1z8XwyhFRKQ/KWGQY7y//TD/89RaAL7+yVkn/eZJn/n4xzmfxjCe3fQKe2sOxCJMERHpZ0oYpEN5VQM/eGgpbUHHhy+axEVnjo5Jv5OGjuPSSefT5oL8dvljGgApIpKElDAIcPSFUlW1TZxxyjCuv7rnT0T0xGdmXkt2ahbrDm3i7T0rYtq3iIjEnxIGAeAPz29g465Khg3J4I7PleD3x/afRk5aNp89/UMAPLTyz9Q118e0fxERiS8lDMKKjYd48rWt+HzGHZ+bc1KDHLtyycTzmFIwkcrGav6w+q9x2YeIiMSHEoZBrqKmkZ/+aTkAn71iKtMnFMRtXz7z8eW5nyPgC/DK9jdYd3Bj3PYlIiKx5bmEwcyuNbOlZrbEzN40s5Ju6m80s0VRn29E1Skys7+Z2dtmtsLMbo/vUSSHtqDjJ39cTnVtM2ecMoyPXzIl7vscnVfEx6Z/AID/XvpHGlub4r5PERE5eYFEBxDJzOYAjwBnOefWm9k1wItmNsM5V3qCZqXOufld9OkDngaed87dbWZ5wAozq3HOPRDrY0gmT7y6mTVbyxmSncZtn52Dv5czOfbVh6ZdwTt7VrCreh+Pr32aL8z+eL/sV0RE+s5rVxi+CbzonFsP4Jx7BjgI3HwSfV4FzAJ+Eu6zGvhv4C4z65/fkB60dU8Vj7y4CYBbP3sm+bnp/bbvgM/PV876PGbGs1teZcvhHf22bxER6RuvJQwLgOiXDiwFLjvJPrc556qi+hwDTD2JfpNWc0sb//XoCoJBx7UXTOTMqcP7PYaJQ8fxwamX4Zzj/vd+T3NbS7/HICIiPeeZhMHMhgJ5QPRUgKXAxC6aZpnZb8NjHhaZ2V1mFvnn8sQT9Nm+bdB55MWN7C49wqjCLD5/1bSExfHJGVdTlDOcvTUHeGTNUwmLQ0REuueZhAHICi+jR8E1AZldtNsE3O+cuxD4OKFbEE9E9dtZn3TWr5ndaGbLzGxZWVlZT2NPGht2VPCXRVvxGdzymTNJT03cMJbUQCpfP+cf8JuP5za/yprSDQmLRUREuualhKEuvIyeBCANOOEsP865zznnloa/lwN3A1eb2ayIfjvrk876dc494Jwrcc6VFBYW9u4IPK6xqZX/enQFzsHHLjmFU8cNTXRITBo6jk+cdg0Av3zvIY401SY4IhER6YxnEgbnXAVQBYyM2jQS2NaLrtrrTg4vt5+gz8i6g8JDz67nQHkd44ty+czl3hm+8eFTr2DqsElUNlTzwLJH9K4JEREP8kzCELYQiJ53oSRcfhwzm2lmN0QVjwovd4eXLwOTzWxIVJ97nHObTi7c5LFhRwXPvLkDv8+49TNnkhLwJzqkDj6fj6+dfT0ZgXTe3buSxTvfSXRIIiISxWsJww+AK8xsGoCZXQUUAb8Mr99rZusiBjUWAP/XzArC29MIPZr5HrA8XOd5YBVwa7hOLnAjcG9/HJAXtLQG+X9PrALg45ecwsRRJ/fK6ngYnj2ML535SQB+u+Ix9tecaNoNERFJBE8lDM655cB1wMNmtgS4E7giYtKmdEIDFdvnT1hDaIDjc2a2CHiT0BMQH3TOtYX7DALXArPN7G1gMfDAYJq06a+LtrK79AhFw7L4xIL4z+bYVxeNP4dzx5bQ2NrET976H5pamxMdkoiIhJnuF59YSUmJW7YselqI5LK/vJav/fg1mluD3HvTuZwxxdsDORtaGvnmyz9g/5GDXDT+HP75rC8wiOfXEhHpbyf8geupKwwSW8457n9iDc2tQS6eM9rzyQJARko6/+fcfyLVn8Line/w2o63Eh2SiIighGFAW7RiL6u2lJGTmcI/XntaosPpsbFDRvFPcz4LwG9WPMbOyr0JjkhERJQwDFC19c385u/rAPiHD84gLzt6Kgpvu2jCOVwy8Txa2lr46VsPUNtc130jERGJGyUMA9SfXt5EdW0zMyYWcOncsYkOp0/+YfYnGT9kNKW1Zdz31m9oC7YlOiQRkUFLCcMAtOfgEZ59YwdmcOOHZybtoMHUQCp3nP9lctOyWXNwAw+vejLRIYmIDFpKGAYY5xz/+7d1tAUdl589zpNzLvRGYVYBt5/3Zfw+P89veY2F295IdEgiIoOSEoYBZtmGg6zYdIis9ACf/0Di3kQZS6cWTuKmkusA+M3yP7H+0OYERyQiMvgoYRhAWlqD/O/fQgMdP335qUk30LEr8yfM45qpC2hzQX7y5gOUHjmU6JBERAYVJQwDyDNvbGd/eR2jCrO5+rwJiQ4n5j53+keYXXQaR5rr+O7iX1DVWJPokEREBg0lDANEdW0Tj74cepfWDR86jZTAwDu1Pp+PW+f9I5Pyx3GwrpwfLPklDS2NiQ5LRGRQGHi/VQapxxZupr6xlTNPHU7JtBGJDidu0lPS+dcL/5kR2YVsr9zNT996gNa21kSHJSIy4ClhGAAOVtTz/Fuhxyivv3p6osOJu7z0XO686GvkpeWwunQD9y/9PUEXTHRYIiIDmhKGAeCRFzfS2ua4aPZoJhQn92OUPTUyu5B/vfBm0gJpvL7rPR5c8Wf0IjURkfhRwpDkdh6o4bXlewj4jeuuPDXR4fSrSUPHccd5NxHwBXhh6yJ+v/ovShpEROJECUOS+/1zG3AOrjxnPCMLshIdTr87feQ0bjvvRvw+P89sWsija/+upEFEJA6UMCSx97cf5r31paSn+vnkZVMSHU7CzCmeya3zbsBnPv664QWeXP9cokMSERlwlDAkKeccDz27HoAPXzSZ/Jz0BEeUWGeNnsXXz/kSZsbj657hifef1ZUGEZEYUsKQpJZvPMSGnRXkZqXykfmTEh2OJ5w7toSbz/piR9LwB41pEBGJGSUMScg5x6MvhSZp+tjFk8lMT0lwRN5x4fizuWXeP+I3H09vWsj/LP+THrkUEYkBJQxJaNXmMjbtriQ3K5UPnDvwpoA+WfPGzOGO879Cij+Fhdte55fvPkRbsC3RYYmIJDUlDEnGOcefwlcXPnzRJDLSAgmOyJvOLD6Nb134VdLD8zT86I1f06hppEVE+kwJQ5JZu62cDTsryMlMGZAvmIqlGcOncPf8b5CTmsXKA+v4t9d+SlVDdaLDEhFJSkoYksyjL20G4EMXTtLYhR44pWAC9y74v4zILmRH5R7uXPgj9lYfSHRYIiJJRwlDElm3rZy128rJSg9wzfkTEx1O0ijKGc69l97OKUPHU1Zfwd2v/Jh1BzclOiwRkaSihCGJPPZy6OrCtRdOIitDVxd6Iy89l3suvpW5o86grqWBexf/nOc2v6rHLkVEeshzCYOZXWtmS81siZm9aWYlXdQdYWbfM7M3zGyRma00s2+aWSCqXlV4e+TnU/E/mtjZuKuCVVvKyEgLcO0FurrQF2mBVG4790Y+dOrlBF2QB1f+mfvf+z3NbS2JDk1ExPM8NcTezOYAjwBnOefWm9k1wItmNsM5V9pJk08AVwDznXNHzGwMsALIAO6JqLfKOTc/zuHH1V9e2wrANedPIDszNcHRJC+fz8d1Z3yECflj+NV7D7No59vsrTnA7efdxNDMIYkOT0TEs7x2heGbwIvOufUAzrlngIPAzSeofwj4sXPuSLj+HuBx4Lp+iLXf7C+r5Z11Bwj4fXxQYxdi4tyxJdx76R0UZg5la8VO7njpu6w68H6iwxIR8SyvJQwLgGVRZUuByzqr7Jx73Dn3aFRxA5AWh9gS5qnF23AOLp4zmvzcwf3OiFganz+G71/+TU4fMY0jTbV8b8n/45E1T9GqSZ5ERI7jmYTBzIYCeUD0M2+lQG/+rJ5H6CpDpJFm9lh4XMRCM/uqmflPItx+U13bxCtLdwPwkfmTExzNwJObls23Lvoqn555LWbGUxte5Nuv/pTyuopEhyYi4imeSRiArPCyKaq8CcjsSQdmtgAYC3wnatNW4C7n3IXADcDXgZ+foI8bzWyZmS0rKyvraexx89ybO2huDTJ3+gjGjMhJdDgDks98fHT6B/i3+beSn5HHpsPbuf3Fe3l953t6ikJEJMxLCUNdeBl9OyENqO+usZmNA34FXOucq4rc5py7xjm3Jfx9J/BD4CvhqxpE1X3AOVfinCspLCzs9UHEUlNLG8+8uQPQ1YX+MH34Kfz48js5s3gm9S0N/OLd3/Ffb/0vNU21iQ5NRCThPJMwOOcqgCpgZNSmkcC2rtqa2XDgb8A/OedW9mB32wCjd7c6+t2rS3dTU9fM5DFDOG1iQaLDGRRy03P4l/O/wpfnfp70QBrv7F3BbS/8B8v2rUl0aCIiCeWZhCFsIRA970JJuLxTZpYPPAt80zm3OFx2Y8T2S83s2qhmo8LL3ScdcZy0BR1PLQ7lSR+9aDJmluCIBg8z45KJ5/KfV9zFtMJTqG6s4Udv3M9/vfW/eheFiAxaXksYfgBcYWbTAMzsKqAI+GV4/V4zW2dm6eH1bOB5QglFmZmVhCd6uimizzHAHWaWGW6TB9wKPOmcO9RPx9Vr771fyv7yOobnZ3Du6UWJDmdQGp49jH+7+Ba+OOvjpAXSeHvPcm55/tss3PYGQRdMdHgiIv3KUxM3OeeWm9l1wMNm1gD4gSsiJm1KJzQAsv3P7W8AZ4c//3qCbl8B5gCvmVkjkAO8Cnw7PkcRG8+8sR2AD14wCb/fa3nd4OEzH1dPvZSzRs/if5c/ysoD63hg2R9ZsvMd/uHMTzE+f0yiQxQR6RemUeAnVlJS4pYti54WIv72HDzCP//oVdJS/Tx4zxVk670RnuCc4+09K/jdysepbqzBzFgw8Xw+NfNactOyEx2eiEgsnPD+t/509aDnwk9GzD9ztJIFDzEzzh07h/s+8G9cPeVSfBgvb3udbzx7D89tfpXWttZEhygiEjdKGDymvrGFV5btAeDq8yYkOBrpTFZqJl+c/XF+fMVdnD5iGnUtDTy48s/c+sJ3eGv3Mo1vEJEBSQmDx7y2fC8NTa3MmFjAhOK8RIcjXRidV8SdF32NO87/MsU5IzhYW8Z9b/+Gb738Q9aUbtCkTyIyoHhq0ONg55zj2fDtiKvP1dWFZGBmzB11BmcWncZrO97mz+ueYXvlbu5d/HOmFU7m4zOu5rThU/VYrIgkPSUMHrJu22H2HDxCfk4a58zUo5TJxO/zs2DS+Vww7iye2/wqf9/0MhvKtvIfi37G1GGT+MSMq5k54lQlDiKStJQweEj71YUrzhlPSkB3i5JRWiCVj0y/kitOuYgXtizimU2vsKl8G/cu/jkT8sfwoVMv5+zRs/H7kuLdZyIiHfRYZRf687HKw9UN/MO9LwPw27suoyAvo1/2K/HV0NLIi1sX88ymhR3vpCjMKuCaKZcyf8I8MlL0unIR8ZQTXgZVwtCF/kwYHnlxI396aRPnnVHMv35hbr/sU/pPc2szi3e+y9ObXqa0NvQW1IyUdOaPn8cVp1xEcc6IBEcoIgJ0kTDoloQHtAUdL78Xeq3FB+aNT2wwEhepgVQum3wBl048j6X7V3fcqnh+y2s8v+U1zhg5ncsmXcCZxTMJ6HaFiHiQEgYPWL25jPKqBkYMzWTmpGGJDkfiyOfzcfbo2Zw9ejY7K/fwwtbFvLHrPVaXrmd16Xry0nK4aMI8Lpl4rq46iIin6JZEF/rrlsQPHl7Km6v387krT+VTl02N+/7EW2qb6li08x1e3f4me2sOdJRPKZjIBePOYt7YOZp6WkT6i8Yw9EV/JAzVtU1c/50XCQYdv7nrcoYN0WDHwco5x5bDO3h1+5u8uWc5Ta1NAPjNx6yiGZw3di5zimdqoKSIxJPGMHjVohV7aW1zzDl1uJKFQc7MmDJsIlOGTeT6Mz/Jsn2reX3Xe6wu3cDy/WtZvn8tKb4As4pmMG/MmZxZPJPMFP2bEZH+oYQhgZxzLAwPdrzs7HEJjka8JD2QxvnjzuL8cWdR1VjDW7uX8c6eFWws38bSfatZum81fp+f04ZPYU7x6ZSMOp1hmUMTHbaIDGC6JdGFeN+S2Ly7ktt+toS87FR+d/cVmqxJulVRX8U7e1fw7t6VbCzfdsz7KsYNGc2skdOZVTSDqcMm6WkLEekL3ZLwovZHKS+eM0bJgvTI0MwhXDXlEq6acgk1jUdYcWAdy/atYXXpenZV7WVX1V7+tvElMgLpzBg+hdNGTGXmiFMZnVukaalF5KQoYUiQxuZWlqzcC8CCs8YmOBpJRrnpOcyfMI/5E+bR3NbCxrKtrDrwPqtK17O35gDL9q9h2f41AOSl5zJj+BSmF05mWuEpjModic+UpIpIzylhSJC31uynvrGVqWPzGTcyN9HhSJJL9adw+shpnD5yGl8AyusqWHtwI2sPbWLdwY0d4yDe2h26xZadmsWpwyaFBlkWTGDi0HGkB9ISexAi4mlKGBLklaV7ALjsbF1dkNgbljWUiyeey8UTz8U5x76aUtaXbWZD2VY2lG2loqHqmCsQPvMxLm8Uk4aOY+LQcUweOo7RecUaByEiHZQwJMDh6gbWbisnJeDj/DNGJTocGeDMjNF5RYzOK+LyyRfhnONQXTmbyrez+fB2tpTvYFf1PnZU7WFH1R7Y/gYAKb4AY/NGMS5/NOOHjGZC/hjG5BXrUU6RQUoJQwIsXrEP52Du9BFkZaQkOhwZZMyMEdmFjMgu5MLxZwPQ2NLI9srdbK/czbaKXWyv2M2B2kNsq9zFtspdx7QfljmUsXnFjMkrZnRuEaNyRzIqd6QmlBIZ4JQwJMDiFaHBjvPPHJ3gSERC0lPSmT58CtOHT+koq29uYGfVXnZW7WFnZWi5r6aU8voKyusrWHFg3TF9FGTkU5w7nKLsERTlDKcoZwQjcwoZnllAwK8fNSLJTv8X97NdpTVs319NVkYKJdP0ciHxrszUDKYPP4Xpw0/pKGsLtlFaW8ae6v3srt7H3ppS9tWUcuDIIQ43VHK4oZK1Bzcd04+ZMSwjnxHZhQzPHsbwrAIKMwsozCpgeFYBQ9Jz8fn0xIaI1ylh6GftVxfOO72YlIAGlEly8fv8HbcgzhlzZkd5W7CNQ3WHOXDkIAeOHAp9ag9SWltOeX0FZeEPhzYd36f5GJoxhILM/I7P0IwhHZ/8jDyGpOeS4tftO5FEUsLQj4JBp9sRMiD5ff7wbYjhx21rbWvlUP1hDtaWUVZ3mEN1FeFlOeX1lVQ31hxNKLqQk5rFkIw8hqTnkJeWS156LnnpOQxJzyU3LZvctBxy03PITcsmzZ+qiapEYsxzCYOZXQvcDTQAfuAbzrkTzs9sZrnAz4EZ4fovA3c651oj6hQBvwaGA2nAI865/4zbQZzAhp0VHKpsYFheOjMmFvT37kUSIuAPUJwzguKczm/BNbe1UNFQRXndYQ7XV1HRUMXhhkoqGqqprK+isrGaqsYajjTXcaS5jj3V3e8zxZ9CTmoWOWnZ5KRmkZ2aRXZqJtlpoWVmSibZqZlkpWaSlZJBZvsyJQO/HiUV6ZSnEgYzmwM8ApzlnFtvZtcAL5rZDOdc6QmaPQjUOefmmlkqsAj4DvCtcJ8+4Gngeefc3WaWB6wwsxrn3ANxPqRjtF9duOjM0fh8+utHBEKTTo3MLmRkduEJ6wSDQWqaa6lqqKaq8QhVjdVUNx6hurGGmqZaqpuOUNN0hJrGWmqaa2kJJyEVDVW9jifNn0pGSjqZKRlkpKSTEUg/dpmSTnog7ZhPWvvSn0p6II3UQCrp/lRSA6mk+lM0q6YMCJ5KGIBvAi8659YDOOeeMbODwM2Erjocw8xOAz4CnB6u32xm9wEPmtn3nHO1wFXALGBBuE61mf03cJeZ/Y/rp7dvtbQGeWP1PiCUMIhIz/l8Poak5zIkvWezoja1NnOkqZaaplpqm+uoba4PL0Pf69o/LaFlfUsDdS0N1Lc00NTWTFNbM1WNNTGLP8WfQpo/lDy0L1P8KaR2fFJJ8QdCZb4UAv4Aqf4UAr6jyxRfgBR/IPS9fekLEPD5CfgC4U/ou9/nJ+Dzh5Z2tMzv8+M3n27XSJ94LWFYAPwwqmwpcBmdJAzh+o3Auqj6GcD5wAvhOtucc1VRdcYAU4GNsQi8Oys3HeJIfQvjRuYwoTivP3YpMmilBVJJCwxlWFbvXvntnKOptYn6lkbqWxuob26gsbWJhtZGGlrCn9ZGmlqbaWhtpLG1icbWJppam2mK/B5OOppam2hua6El/PEKv/nCyYP/mETCbz587d87yvz4fD587dvNhz+8bhFlvqjvPrOI76F1CycrR7cbxtFtoXUL920d65FtjNA2g2PaWOTSCNfj+G0AhPo+fklHMtXef/tae1/Hfqej76Pf21t3tn60XXRb2mM9uvG49aM9RfRjxrDM/H65iuWZhMHMhgJ5wIGoTaXAB07QbCJwMOoqQWnEtvZlZ322b+uXhGFRxO0IEfEmMyM9JZ30lHSGMiQmfQZdkJa2VprDSURzWwvNrS00tzWHPy20BFtpaWuhOVyvNdhKc1tLeNlKa3udYCutba20BttoCbbQGmyjNaKs1UWsuzbagqFP+7a2YBtBF6TNBWlrCwLeSWKk7x786E/7ZQZWzyQMQFZ42RRV3gRkdtGms/pEtOlJnQ5mdiNwI8DYsbF5z4Nzjpq6JszgotlKGEQGE5/5wlc8UslJdDCEEphgMEira+tYHpNIBI9dD4bLOlsPfRxtLrwebP/uIraHPs45HK6jTXtZ0AWPKe+oFwziwvFG1nPOEeRoPeeOlhG5rX37cetAZFsIL4O0/+kZrkXQudCai2wX3n7Mdzr2RUQf7W1C60fbHbMeXomILmI9rKPN0baR5ZFXLeLJSwlDXXgZ/cq8NKC+izad1SeiTR0QfeMzuk6H8EDIBwBKSkpiMr7BzLj3y+dxuLqBgjzNwy8iieMzHz6/j4CnfvxLMvDM0F3nXAVQBYyM2jQS2HaCZtuB4XbsCJ729tsi6nTWZ2SdfqFkQUREkpVnEoawhUBJVFlJuLwzLxMa4Dgjqn4D8GZEnclmNiSqzh7n3PHTzomIiMhxvJYw/AC4wsymAZjZVUAR8Mvw+r1mts7M0gGcc+8DfwXuCG9PAb4B3Bd+pBLgeWAVcGu4Ti6hMQr39tMxiYiIJD1P3cRyzi03s+uAh82sfabHKyImbUonNFAx8hbE9cAvzGxpuP5C4J6IPoPh2SN/bWZvh/t4oL8nbRIREUlm1k/zFiWlkpISt2zZCWelFhERGWhO+MiF125JiIiIiAcpYRAREZFuKWEQERGRbilhEBERkW4pYRAREZFuKWEQERGRbilhEBERkW4pYRAREZFuaeKmLphZGbArhl0OA8pj2F8i6Vi8Z6AcB+hYvGqgHMtAOQ6I/bGUO+eu7GyDEoZ+ZGbLnHPRL9dKSjoW7xkoxwE6Fq8aKMcyUI4D+vdYdEtCREREuqWEQURERLqlhKF/DaQ3ZOpYvGegHAfoWLxqoBzLQDkO6Mdj0RgGERER6ZauMIiIiEi3AokOYCAxs2uBu4EGwA98wzm3rIv6ucDPgRnh+i8DdzrnWvsh3C714Vg2AqVRxX91zv0sflF2z8xSgW8DdwCTnXM7u6nv5XPS22Px3Dkxs2uAfwbSgCygEvgX59yabtp57rycxLF48bxcCNwC5BP675sP/G9XMXn0nPTlODx3PqKZ2dcI/be+2Dm3qIt68T0nzjl9YvAB5gC1wPTw+jXAYWBkF23+Avw+/D0VeAv4XpIey6JEx91JTOOBt4GHAAeM70Ebr56TvhyLF89JOfDZiPUfAGXAiGQ7LydxLF48L78G7olYnwW0Adck2Tnpy3F47nxExVdMaD4gB8zvpm5cz0nC/2MMlA/wBPBkVNl64D9OUP+08D+AmRFlnwTqgexkOpbw9kWJPgcn+G88GZjfk1+yHj8nvToWD5+Tv0StF4aP5/PJdl76ciwePi/TgZyossPALUl2Tnp1HF49H1HxPQl8ubuEoT/OicYwxM4CIPqS/VLgsi7qNwLroupnAOfHPLre6e2xeJJzbp1zbmsvmnj2nPThWDzJOffRqKKG8DKti2aePC99PBZPcs6td84dATAzn5n9E9AE/PkETbx6Tnp7HJ5mZh8EWoAXelA97udEYxhiwMyGAnnAgahNpcAHTtBsInDQhdPAiPrt2xKij8cCkGVmvyX0V3AQWAj8p3OuMS6Bxocnz8lJSIZzMo/QD7m/d1EnWc5LT44FPHxezOwu4KtABXCVc27fCap6+pz04jjAo+fDzLKA7wJX0LMkNO7nRFcYYiMrvGyKKm8CMrto01l9umjTH/pyLACbgPudcxcCHweuInRrI5l49Zz0lafPiZkZoYG1dznnDnVR1fPnpRfHAh4+L865e4EiQr+oFpvZuSeo6ulz0ovjAO+ej/8Afu2ci/7j7UTifk6UMMRGXXgZnQWmEbp/dKI2ndWnizb9oS/HgnPuc865peHv5YR+eF5tZrPiEWScePWc9EkSnJPvAbuccz/ppl4ynJeeHovnz4sL+SOwhNBAzs54/pz08Dg8eT7MbDZwNqFBnD0V93OihCEGnHMVQBUwMmrTSGDbCZptB4aH/zKJrE8XbeKuj8fSmfa6k2MQVn/x5DmJIc+cEzO7BZgGfKkH1T19Xnp5LJ1J+HkJP7IbbT2hx/M648lz0ofj6EzCzwehJ9MygFfNbBHwaLj8PjNbZGZTO2kT93OihCF2FgLRbwwrCZd35mVC/yAi/yGXEBo49WbMo+udXh2Lmc00sxuiikeFl7tjHFs8efmc9IqXz0k4rquATznnWs1sopkt6KKJZ89Lb4/Fw+dledQvGgg9zneie/9ePSe9Og6vng/n3H845850zs13zs0HPh3edEu4bFMnzeJ/ThL9yMhA+RCau+AIMC28fhWhATcjw+v3Ehq9mh7R5i/AQ+HvKcAbeOOZ/14dC6FH/TYDBeH1NOAZ4F3A74HjmU8njyIm0znp7bF49ZwQ+sG3A7iI0A+zEuAm4N+T7bz05Vg8fF52Al+PWJ9D6BfNN5LsnPTqOLx6Pjo5rvFEPVaZiHOipyRixDm33MyuAx42s/bZEa9wzrWPUk0nNPAkMvu9HviFmS0N118I3NN/UXeuD8eyhtAgoefC9bOBVcA/OOfa+jX4COHLky8BQ8JFj5rZfnf0cbikOSd9OBZPnhPg94SezloUVf7tiO/Jcl76cixePS/fAm4ws88QmugoA7gNuD+8PVnOSW+Pw6vno4OZ3QecE169z8y2OOc+QQLOiV4+JSIiIt3SGAYRERHplhIGERER6ZYSBhEREemWEgYRERHplhIGERER6ZYSBhEREemWEgYRERHplhIGERER6ZYSBhHxFDP7npn9PNFxiMixlDCISL8ws/80sxd6UPUq4Ll4xyMivaOEQUT6y1zgva4qmFkxcArHv59BRBJM75IQkbgysxSgjtDb89qtd87N6KTuDcCHnXPXdNHfZ4FbgVpCL+B5zTn3rdhGLSLRlDCISFyZmQ+YDSwDzgZ2A03OucpO6v4FeNk5d3/0tvD24nD7Kc657WZWCGx0zhXE7QBEBNAtCRGJM+dcECgCjgBLnXOlJ0gWUoAFdD1+YQSh1/aOD/ddRmjMg4jEmRIGEekPs4HVrutLmhcCe5xzu7qoswr4PbDQzF41sxuBdbELU0RORAmDiPSHWcDKbup0+3SEC/kCMBNYDnwXWGVmQ2IQo4h0QQmDiPSHM4A13dS5Cni2qwpmNsrM5jnn3nfO3QHMAIqBS2MTpoiciBIGEekPAeBUMyvu7GqAmU0g9Iv/zW76OQX4oZkFwus+wICtMYxVRDqhhEFE+sOdwKeBvcD3O9l+NaGnI1q66WcjsA1428wWAU8DNzvnVscwVhHphB6rFJGEM7PngCecc79NdCwi0jldYRARL1hE6GqBiHiUrjCIiIhIt3SFQURERLqlhEFERES6pYRBREREuqWEQURERLqlhEFERES6pYRBREREuqWEQURERLr1/wFG9IqTBcn2qQAAAABJRU5ErkJggg==\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 generalized midpoint rule. It can be considered as a weighted average between an explicit forward Euler and an implicit backward Euler scheme. For $\\alpha > 0$ the update proceeds implicitly with the *recursive* update formula\n",
    "\n",
    "$y_{n+1}^{i+1} = y_n + [(1-\\alpha) f(y_n) + \\alpha f(y_{n+1}^i)] (t_{n+1} - t_n)$\n",
    "\n",
    "We define the corresponding Euler update as\n",
    "\n",
    "$\\Delta y_\\text{MP}^{i} = y_{n+1}^{i+1} - y_n = [(1-\\alpha) f(y_n) + \\alpha f(y_{n+1}^i)] \\Delta t$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Mid_Point_Update(y_i,y_n,al,dt):\n",
    "    return ((1.-al)*dy_dt(y_n) + al*dy_dt(y_i))*dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, the first step is a forward Euler step. After that, a mixture of BE and FE schemes occurs.\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 generalized Midpoint 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 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,alpha):\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]+Mid_Point_Update(y_old,y[-1],alpha,dt) #one backward Euler step\n",
    "            i += 1\n",
    "            #if change is below tolerance, above maximum iterations, or scheme is fully explicit, stop iterations\n",
    "            if (np.abs(y_new - y_old) < absolute_tolerance or i > max_iter or alpha == 0.0): \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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d50fd1a61718451fae512da9580c877a",
       "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",
    "          alpha=widgets.RadioButtons(options=[0,0.5,1], value=1, description='$\\\\alpha$'))\n",
    "\n",
    "def plot_CN(dt=0.1,alpha=1):\n",
    "    time_n, y_n, iter = integrate(dt,alpha)\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='numerical',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": [
    "We can also plot the error:"
   ]
  },
  {
   "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$ and $\\Delta t = 0.3$?\n",
    "- Can you recover the forward and Backward Euler solutions?\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.4"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}