03b_Generalized_Midpoint_Picard_Nonlin.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 3 - Generalized midpoint rule, Picard iterations\n",
    "\n",
    "The ODE $\\dot{y} = k (y_\\text{eq} - y)^2$ 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": "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):#assuming y0=0\n",
    "    return k*t*y_eq**2./(1.+k*t*y_eq)"
   ]
  },
  {
   "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)**2"
   ]
  },
  {
   "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,5]$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "time = np.linspace(0,5,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": "iVBORw0KGgoAAAANSUhEUgAAAgwAAAF9CAYAAABoEDD2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABK+0lEQVR4nO3deZydZX3//9fn7LOvSSbbJCSBhIRAICEQBQwIgoi4W6tWaUvh6w+rWLGVCnyrRYVWW6pfSou1oFWKKIqCCLIYAsiSlayErGSdyWQmk9nXc/3+OGcmJyeTzCQ5c+77zLyfj8d53Ofc93Xf92eG5bznvq77us05h4iIiMjxBLwuQERERPxPgUFEREQGpcAgIiIig1JgEBERkUEpMIiIiMigFBhERERkUCGvC/Czq666yj311FNelyEiIpItdqwNusJwHAcOHPC6BBEREV9QYBAREZFBKTCIiIjIoBQYREREZFAKDCIiIjIoBQYREREZlG6rFBGRUSsej3PgwAEaGxvp7e31upxhF4vFmDRpEuFw+IT3VWAQEZFRa/fu3ZgZU6dOJRwOY3bMaQhynnOO+vp6du/ezWmnnXbC+6tLQkRERq3W1lYmTpxIJBIZ0WEBwMyoqKigo6PjpPb3TWAws2vM7Ekze87MXjWz35nZ2UPYr9jMHjSzZWa20szuNrNQWpvxZvZrM3sl2eaW4ftJREQklwQCvvkqHHanEor89Ft6EPiJc+7dzrkLgTeA58xs3BD2CzrnzgcuBC4GvtG30cwCwOPAGufcIuBS4HNmdkPmfwQREZGRyU+BYalz7qGUz98FKoH3HGsHMzsL+BDwTwDOuS7gHuBmMytMNrsamJc8Hs65Q8B/ArfZSL/+JCIikiG+CQzOuQ+nrWpPLqPH2e1yoANYl7JuGZAHXJTSZqtzrjGtzWRg5snWe6L+uHMFP171C3Yc3J2tU4qIiGSMbwLDABaRCAO/OU6baUCtc86lrKtJ2da33Je2X3qbfmZ2g5ktN7PldXV1J171MazYu4Yn3nqOtxsVGEREJPf4MjAkuwpuB25zzu0/TtMCoDNtXd/n/BNo0885d79zboFzbsGYMWNOrPDjKIomekiau1oydkwREZFs8WVgAL4FvO2c++4g7Vo5usui73PbCbQZdsXJwNDUqcAgIiKnbtKkSfzLv/zLEevWrl1LLBZjw4YNGT+f7yZuMrObgdnAR4bQfBsw1swspVuiKrncmtLmqrT90tsMu6KIAoOISC54/5d/7cl5H//uB06o/aJFi1i2bNkR626++Wauv/56Zs+encnSAJ9dYTCz60nc1fBx51yPmU0zs8uPs8szJAY4zklZt4DEgMmXU9rMMLPStDa7nHObMlb8IIqiBQA0KzCIiEgGpAeGxx57jFWrVvH1r399WM7nmysMZvYJ4GvAdcDc5B2P84HxwLPJNncCHwQWOOc6nHPrzexXwFeAz5pZGPgicI9zru+b+XfAauBLwP81s2LgBuDO7PxkCcXRIkCBQUTE7070L32vXHjhhXz5y1+moaGBgoICbrnlFu644w4qKiqG5Xy+CQzA/5CoZ0na+tSoFCMxUDF1/oTrgO+b2TIgSCJc3NG30TkXN7Nrgf8ws1eSx7jfOXd/pn+A4+kbw9Dc2ZrN04qIyAg1f/58IpEIy5cvZ9WqVYRCIW666aZhO59vAoNzbtBHZznnbgFuSVvXBHx2kP32AteeUoGnqK9Lokl3SYiISAZEo1HOPfdcHn/8cX70ox/x0EMPndRTKIfKV2MYRrLCSCIwtHS2Eo/HPa5GRERGgkWLFnHvvfeyaNEirrnmmmE9lwJDlgQDQQoi+TgcLd1Zu5tTRERGsHnz5hEIBI66vXI4KDBk0eFxDOqWEBGRU/eTn/yEG2+8kTlz5gze+BT5ZgzDaFAcKWQf+xUYRETkpMXjcerq6njwwQdZt24djzzySFbOq8CQRUWa7VFERE7R0qVLueyyy5g5cyaPPvooZWVlWTmvAkMWKTCIiMipWrx4sSeD5zWGIYuKNIZBRERylAJDFmnQo4iI5CoFhizqf2KlJm8SEZEco8CQReqSEBGRXKXAkEVFydkeNehRRERyjQJDFmkMg4iI5CoFhiw6/IhrPbFSRERyiwJDFuWFYwQtQHtPB9293V6XIyIiMmQKDFlkZikDH3WVQUREcocCQ5ZptkcREcmknTt3UlhYyKpVq4b1PJoaOsv6Bz5qLgYREcmACRMmsHr1aiZPnjys51FgyDLNxSAiIpkUCoWYMWPGsJ9HXRJZVhxRl4SIiOQeBYYs0xgGERHJRQoMWVYUTcz2qC4JERE5FZMmTeJf/uVfjli3du1aYrEYGzZsyPj5FBiyTLM9iohIJixatIhly5Ydse7mm2/m+uuvZ/bs2Rk/nwY9Zln/bI+6S0JExJc+/rPPeXLeR/7kvhNqv2jRIv793/+9//Njjz3GqlWreOSRRzJdGqArDFnXP4ahQ4FBRERO3oUXXsjWrVtpaGigs7OTW265hTvuuIOKiophOZ+uMGRZ3xiGJl1hEBHxpRP9S98r8+fPJxKJsHz5clatWkUoFOKmm24atvP5KjCYWQT4BnALMMM5t+M4ba8DvgrUpG1aBHzCOferZLtGYHVam/uccz/LSNEnqChyeGpo5xxm5kUZIiKS46LRKOeeey6PP/44P/rRj3jooYcIh8PDdj7fBAYzmwr8L/AWEBzibnc55x5MOcY0YAXwVEqb1c65xZmp8tRFQxGiwQidvV109HSSF455XZKIiOSoRYsW8W//9m9cccUVXHPNNcN6Lj+NYSgE/gx4YIjtH0u+Un0WeNg51565sjJPd0qIiEgmzJs3j0AgcNTtlcPBN4HBObfOObflBNo3Ouca+z5b4tr+Z4D/HobyMkqTN4mISCb85Cc/4cYbb2TOnDnDfi7fdElkwKVAi3NuWdr6KjP7GTAe6CJxVeI+51xvluvrp8AgIiInKx6PU1dXx4MPPsi6deuG7TbKdCMpMFwH/HCA9VuA25xzm5PjJH4PnAkMOJTUzG4AbgCorq4elkL1ACoRETlZS5cu5bLLLmPmzJk8+uijlJWVZeW8IyIwmFkR8H7gb9K3OeeuSXm/w8zuBn5gZrc75xoGaH8/cD/AggUL3HDUWxxJTg+tWytFROQELV68mHg8nvXz+mYMwyn6OPCcc+7AENpuBQyYNrwlHVtxLDHbo7okREQkV4yUwHAdA3RHmNm7zezatNUTk8udw13UsRTpEdciIpJjciYwmNmdZrbOzGJp62cAU4GnB9htMvAVM8tPti0BvgQ86pzbP8wlH5OeWCkiIrnGN2MYkrM8/h4oTa562Mz2Ouc+nPwcA/JJdCekug74sXNuoA6d54DzgD+YWQdQBDwPfD2z1Z8YzcMgIiK5xjeBwTnXBSw+zvZbSEwZnb7+tuPsswv4Qibqy6T+J1Z2tnpciYiIjKZp+p07+bH8OdMlMZLoAVQiIv4QDodpb/f15MAZ1d3dTSh0ctcKFBg8UJi8rbKls9WTW2NERCRh7Nix7Nmzh7a2tlP66zsXxONxamtrKSkpOan9fdMlMZoEA0EKIvm0drXR0t3WP6ZBRESyq7i4GIC9e/fS3d3tcTXDr6CggMrKypPaV4HBI8WRQlq72mjubFFgEBHxUHFxcX9wkGNTl4RHdKeEiIjkEgUGj+gBVCIikksUGDyiwCAiIrlEgcEjxf2BodnjSkRERAanwOCRsrzEbS0N7Y3eFiIiIjIECgweKc8rBaCh/ZC3hYiIiAyBAoNH+gLDwbZGT+sQEREZCgUGj/RfYeho9LQOERGRoVBg8EhpLDFJSGNHE73xXo+rEREROT4FBo+EgiFKokU45zikOyVERMTnFBg81N8toXEMIiLicwoMHtKtlSIikisUGDzUf6eEbq0UERGfU2DwUHl+KaArDCIi4n8KDB4qi6lLQkREcoMCg4f6rjCoS0JERPxOgcFDZbFSQFcYRETE/xQYPKQxDCIikisUGDxUFCkgFAjR1t1OR0+n1+WIiIgckwKDh8ysfy4GjWMQERE/U2Dw2OHHXDd6WoeIiMjx+CowmFnEzO4ysx4zmzqE9m+a2ZK01xfT2hSb2YNmtszMVprZ3WYWGrYf4gQdnryp0dM6REREjsc3X5zJgPC/wFtAcIi71TjnFg/S5kGg1Tl3vplFgCXAN4C/P6lCM0zTQ4uISC7w0xWGQuDPgAcydUAzOwv4EPBPAM65LuAe4GYzK8zUeU6FHkAlIiK5wDeBwTm3zjm3JcOHvRzoANalrFsG5AEXZfhcJ6U/MHRo0KOIiPiXb7okTlKBmf03MAOIA88C33HOdSS3TwNqnXMuZZ+alG2eK++7S0JXGERExMd8c4XhJG0C7nPOXQJ8FLga+EXK9gIgfYKDvs/5Ax3QzG4ws+Vmtryuri7T9R5Fd0mIiEguyOnA4Jz7tHNuWfL9AeB24H1mNi/ZpBWIpu3W97ntGMe83zm3wDm3YMyYMcNQ9ZHKUrokjrwQIiIi4h85HRgGsDW5nJFcbgPGmpmltKlKa+upaChCQTiP3ngvzZ0tXpcjIiIyoJwNDGY218yuT1s9MbncmVw+Q2KA45yUNguAduDl4a1w6NQtISIifpczgcHM7jSzdWYWS66qAP7WzCqS26PArcBrwAoA59x64FfAV5JtwsAXgXucc775c76/W0LTQ4uIiE/55i6J5KRKvwdKk6seNrO9zrkPJz/HSAxU7OteWENigOOTZtZOYh6H1cBfOOd6Uw59HfB9M1tGYkKoZ4E7hu8nOXG6wiAiIn7nm8CQnFRp8XG23wLckvK5gSHM1uicawI+m4ESh015vmZ7FBERf8uZLomRrCxWCuiJlSIi4l8KDD5Qnl8K6AqDiIj4lwKDD+iJlSIi4ncKDD6gJ1aKiIjfKTD4QGm0mIAFaOpsobu32+tyREREjqLA4AOBQIDSWDEAjR1NHlcjIiJyNAUGn1C3hIiI+JkCg09o8iYREfEzBQafOHynhOZiEBER/1Fg8Im+Lon6toMeVyIiInI0BQafGFtQAUBda4PHlYiIiBxNgcEnxhZUAlDbWudxJSIiIkdTYPCJcYXJwNByAOecx9WIiIgcSYHBJ4qjRcRCUdq622ntavO6HBERkSMoMPiEmTGuv1vigMfViIiIHEmBwUfG9ndLaByDiIj4iwKDj4wrHAMkxjGIiIj4iQKDj6hLQkRE/EqBwUfGqUtCRER8SoHBR/q6JParS0JERHxGgcFHxuSXYxgH2g/S09vjdTkiIiL9FBh8JBQMUZFfhnOOA22aIlpERPxDgcFn+sYx1KhbQkREfESBwWf67pTYr2dKiIiIjygw+MzYlGdKiIiI+EXI6wJSmVkE+AZwCzDDObfjOG3HAV8ELgF6gBLgEeCfnXM9Ke0agdVpu9/nnPtZJmvPlCpN3iQiIj7km8BgZlOB/wXeAoJD2OVjwJXAYudcs5lNBlYCecAdKe1WO+cWZ7ba4TNWkzeJiIgP+alLohD4M+CBIbbfT+JqQjOAc24XiSsMnxqe8rKjb9Djfj3mWkREfMQ3Vxicc+sAzGzSENs/MsDqdiCaybqyrTBSQH44j7budpo7WyiOFXldkoiIiK+uMGTCIhJXGVJVmdnPzGypmT1rZp83s6F0eXhCj7kWERE/GjGBwcwuB6pJDJpMtQW4zTl3CXA98AXge8c5zg1mttzMltfVeXNro+6UEBERvxkRgcHMpgD3Adc65xpTtznnrnHObU6+3wHcDXzOzMoHOpZz7n7n3ALn3IIxY8YMb+HHcPgx15qLQURE/CHnA0Py9spfA9c751YNYZetgAHThrWwU6AuCRER8ZucDgxmVgY8AdzqnHshue6GlO3vNrNr03abmFzuzE6VJy71TgkRERE/yJnAYGZ3mtk6M4slPxcCvwOeBerMbIGZLQBuTNltMvAVM8tP7lMCfAl41Dm3P7s/wdD1BQZdYRAREb/wzW2VyVkefw+UJlc9bGZ7nXMfTn6OAfkkuhMgMcvjBcnXV49x2OeA84A/mFkHUAQ8D3w94z9ABlXklxOwAA1tjXT3dhMOhr0uSURERjnfBAbnXBew+DjbbyExZXTf528C3xzkmLtI3BWRU0KBIJX5ZexvraeutZ4JxVVelyQiIqNcznRJjDb9d0qoW0JERHzglAKDmV2QqULkSP13Smjgo4iI+MCpXmH4eUaqkKNo8iYREfGTQccwmNlAz2yAxODDASc/klPX95jrfc21HlciIiIytEGPl5N4imRL2noDLsl4RQLApJLxAOxq2udxJSIiIkMLDEuAZufc0vQNZrYm4xUJAFWFYwkGgtS11tPR3UEsHPO6JBERGcUGHcPgnPvwQGEhue2KzJckkLi1ckLROAB2N9V4XI2IiIx2Jzzo0cw0KUCWTC6ZAMCuQ3s9rkREREa7k7lL4vcZr0IGNLlY4xhERMQfTiYw2OBNJBP6rjDs1hUGERHx2MkEBpfxKmRAh7skdIVBRES8pamhfWxcQSXhYJj69oO0dbV7XY6IiIxiCgw+FggEmNh/p4SuMoiIiHdOJjD0ZrwKOSbdKSEiIn5wwoHBOXfucBQiA1NgEBERP8hIl4SZFZvZh8zsrEwcTw6bpFsrRUTEB04qMJjZt8yszszON7N8YBnwP8CrZvaZjFY4ylX331qpwCAiIt452SsMi4EznXPLgE8BZcBUYAZwU0YqEwAqC8qJBiMc7DhES2er1+WIiMgodbKBod05dyD5/hPAA865A865GqAtM6UJQMACKd0SGscgIiLeONnAUGRmU8zsUuBdwIMAZhYE8jNUmyRpAicREfHaUB5vPZB7gC0kAsf/OOc2mtmFwN3AugzVJkmTSpJXGHSnhIiIeOSkAoNz7iEz+wMwzjm3Orl6J3AH8GaGapOkycnAoMmbRETEKyd7hQHn3D5gX8rnvYD+BB4GmotBRES8pqmhc0BFXhl54RhNnS0c6mjyuhwRERmFFBhygJkdvlNCAx9FRMQDvgsMZhYxs7vMrMfMpg6hfbGZPWhmy8xspZndbWahtDbjzezXZvZKss0tw/YDDJO+bgmNYxARES8MOTCY2T1mdoWZRYarmGRAeAEYDwSHuNuDQNA5dz5wIXAx8I2UYwaAx4E1zrlFwKXA58zshsxVPvwmJ68w7NQ4BhER8cCJXGFoAf4JOGBmvzGzG81scobrKQT+DHhgKI2Tz674ULIunHNdJG75vNnMCpPNrgbmAd9NtjkE/Cdwm5lZBmsfVn1XGHY27vG4EhERGY2GHBicc7cln1Q5E/g18B5gnZmtTXYDXJL8a/6kOefWOee2nMAulwMdHDn3wzIgD7gopc1W51xjWpvJJH6WnDCtrBqAHY276InrCeMiIpJdJ/N4633OuR865z4CVAJfJNF98B9AvZk9YmbXmVlFhmsdyDSg1jnnUtbVpGzrW6Z3/Ke38b3CaAHjCsfQ1dvNbnVLiIhIlp3qFYFu59zzzrlbnHOzgXNJjEH4GPDnmShwEAVAZ9q6vs/5J9Cmn5ndYGbLzWx5XV1dxgrNhBnlUwDY0vC2x5WIiMhok9G7JJxzO5xz9zrn3uec+04mj30MrUA0bV3f57YTaNPPOXe/c26Bc27BmDFjMlZoJswonwrAloYdntYhIiKjj+9uqzxB24CxaYMXq5LLrSltqjhSepucMKNiKgBb63d4WoeIiIw+uR4YniExwHFOyroFQDvwckqbGWZWmtZml3NuUzaKzJSppZMJWICdTXvp6EnvZRERERk+ORUYzOxOM1tnZjEA59x64FfAV5LbwyQGYd7jnGtJ7vY7YDXwpWSbYuAG4M7sVn/qoqEI1SUTcM6x/eBOr8sREZFR5KQCg5nNynQhyeNGzGwJibkUAB42s1+mNImRGKiY2gVxXXLfZcBrwB9JPDUTAOdcHLgWONfMXiExKPN+59z9w/EzDLe+cQxbNfBRRESy6GSfVrnKzO4H/sE5dzBTxSQnXlp8nO23ALekrWsCPjvIcfeSCA05b0bFVJ7d9hJbNI5BRESy6GS7JBaSGDew2cz+2syGOo2znKLp/bdW7vC2EBERGVUGDQxmdk5ybEA/59xa59zlwPXAF4C1ZvbeYapRUkwqHk80GGF/az1NHc1elyMiIqPEUK4wvERiGuWjOOceI3Gl4Uckxhs8OVzjGyQhGAgyrTwxTfTWgxrHICIi2TGUwDAbON43Uz6wgkRouBJYY2bfM7OSDNQnA5jeN4GTxjGIiEiWDDro0Tm3K32dmd0MnJ98TQe6SNy6+G/J5aeBDWb2Yefca5krV0BTRIuISPad7F0SXwZeAe4DXgVWJO9w6PNjM/s74L85clIlyYDUKaKdc+TQU7pFREaN3t44nd29dPccXnZ19ybed/et66WrO053Ty+dyWVXd5yunkSbrp7Efn379i27ehJtu3vi3PrZhYyvLBj2n+ekAoNzbsAxDWkeAL51MseX4xtTUEFRtJDmzhbqWusZW1jpdUkiIr4Xjzu6enrp7Ep8afctuwZ4n/hij/e37/uiP7yMH7Gu70u8q/vwvvG4G7yoDGjv7MnKeU72CsNQ1AGXDePxRy0zY0b5VFbtW8eWhrcVGERkROiNOzq7eujo6qWjM7lMfu7s6j28rauXzu6e5Lrew9u7U9p2DxwMsilgEAkHCYeCRMMBwuEg0XCQcChAJBwk0rcMB4mEA4RDiWUklNgWHqBtOJTYHg4H+t9PGDP8VxdgGAODc86RmFVRhsH08inJwLCDd1TP97ocERlFnHN098Rp7+yhPfnF3t7RQ3tX8nPy1Z784m/v6qGj88gQ0J4aCDoTy+6e+LDXHgknvryjyS/qWCREJBwgGgkmtwWPfN//hZ7cL5LyOXT4yz5yRNvE52DARlSX8XBeYZBh1D+OoX67t4WISE5wztHVE6etvZu2zh7aOrpp7+yhrSPxau9IrG/v7KG9o+eI9+2dKZ+TYaB3GC63m0E0HCQWDRGLJL7Mo5EgsUiQaDixLhpJbI+Gk+uTbfo+x1I+R5PtU7/4A4GR8wWebQoMOer0vkddN7xNV283kWD4+DuISM5yztHR1UtbRzct7d20tffQ2tFNa3v34XUdiRDQ1pHY1ve5tS8MdGT2Sz4UNPKiIfKiIWLJZV4kRF7s8Jf94W2JL/a85Jd9XiRENBpMLCNB8qKHv+RH0l/kI40CQ44qihYypXQSbzfu5q0D2zhr3EyvSxKR4+iNO1rbu2lp66KlvZuWtm5a2lPfJwJAS3sXrcn3re09ifUd3RkZQBcOBSiIhcmLhsjPS3yh50eTn2OJV16sb32IvGiYvFjf+5RgEE30y8voktHAYGYlzrlDmTymHNtZY2fyduNu1u3fpMAgkiW9vXFa2rtpau2iua2L5tYumtu6E+/bumhJvm9p66a5PbGtta2L1o5TG8keCQcpzAuRHwtTEAtTkJd45cdCFMTC5Ocll7GUdbFQ/+f8WEhf8nJKMn2F4aCZbQFudc49muFjS5qzxs3kt289x/raTTDX62pEco9zjtaOHppaO2lq6eJQSydNrV39r0Otic/NKeta2rtP6lxmUBALU5gfpjAvTGFehIL+92EK8yMU5B3+3Pc+PxkOwqGTfVagSGZkOjC8BzgTuNvMgs65RzJ8fElx5pgZBCzAloYdtHd3kBeOeV2SiOe6untpbO6ksaWzf3mopZNDyUCQ+rmptZOe3hO71N/3xV9UEKE4P0JRQYSi/DBFfe/zEtsK8yKJcJDclh8LE9SAO8lhJxUYzGyWc+7N9PXOuWeBZ83sJxzjgVWSOfnhPKaXVbO5YQdvHtjCuePP8rokkWHRG3ccaumkoamDxubE8mBzB41NnRxMhoKDTR00tnTSdoKX/vOiQYoLopQURiguiFJcEKGkMLFMfRXlJ5aF+RF98cuodLJXGFaZ2f3APzjnDqZvTK47ar1k3pxxM9ncsIN1tZsUGCTnxOOOptYuGpo6qD/Unlx2HLE82NTBoZZOhjrmLxQ0SgqjlBZFKU1ZlhQmQkFJ3/tkSIiE1a8vMhQnGxgWAv8KbDazrwP/7pzL7hRaAiQGPj628WnW73/L61JEjhCPOxpbOjnQ2E5dYzsHkq/6Q4lwcOBQBw2HOujpHdpkPSWFEcqKYpQVRSkrPnJZWhSlrChGaVGUwrywbs0TGQYn+yyJtcDlZvZB4J+Bz5nZl51zv8tkcTK4mZXTCQaCbD+4i5auVgoj2ZkiVKSzu5f9DW3UHWxn/8E26hrbqetftlN/qH1I4wMK88JUlMSoKMmjvDhGRUmM8pIY5cWJV18Q0KA/EW+d0qBH59xjZvYk8CXgYTN7GfibgcY3yPCIhiKcUTGNjXWb2bB/MwsnzfO6JBkhunvi1DW2UVvfRm1DGzX1rew/2M7+hjZqD7bR2Nw56DGKCyJUluYxpjSPytI8KkpiVJbmUVmSR0VpIhDEIpoORiQXZOK/1HxgBfAj4CZgjZn9B3C75mTIjrPGnsHGus2s279JgUFOSHtnDzX1rew90Mq+A63U1B9eHmhsP+64gVDQqCzNY2xZPmPK8hhTms/YsrzE+7J8KkvziGp8gMiIcbJ3SdwMnJ98TQe6gNXAvyWXnwY2mNmHnXOvZaJQObazxs3k5+t/m5iPQSRNd08vew+0sreuhT11iWXf54PHuUpgBpWleVRV5DOuPJ9x5QXJZeJVVhzT3QIio8jJXmH4MvAKcB/wKrDCOdeVsv3HZvZ3wH8Dc06tRBnMjPKpRIJhdjXto7GjidJYsdcliQeaWrvYVdvM7v3N7N7fknw1s7+h7ZhXCkLBAOMr8xlfUUhVZT4TKgoYV1HA+MoCxpblaWZAEel3soMehzLHwgPAt07m+HJiwsEwsypnsKZ2Ixv2v8U7qhd4XZIMo+a2Lnbsa2JnTTM7a5rYVdvCztomDrV0Ddg+YDC+soCJYwqZMCaxnFhZyPjKAipK83SVQESGZDhHG9UBlw3j8SXFWeNmsqZ2I+tqNykwjBDdPb3s3t/C9r2H2L63ibf3NfF2TRMNTQN3I+RFg0wcW8TksYVMHlfEpLGFTBpbRFVFge4wEJFTNmyBwTnngBdOdD8zuxa4HWgHgsAXnXPLj9H2OuCrQE3apkXAJ5xzv0q2ayQxtiLVfc65n51ofX511tjEw6fW7tc4hlzU0tbFtr2H2LaniW17Gtm+t4ldtc0DPo44GgkypaqI6nHFVFcVJV7jiqksjWn+AREZNr66n8nM5gMPAQudcxvM7BrgaTOb45xLDwV97nLOPZhyjGkk7tp4KqXNaufc4mEq2xdOK5tMQSSf2pY69jTVMLG4yuuS5Bha2rrYsruRLbsPsWVXI1t2N1Lb0HZUOzOYOKaAqeNLOG1CMVPGFzN1fDFjy/IJqBtBRLLMV4EBuBV42jm3AcA594SZ1ZK4XfP2Ado/NsC6zwIPO+fah61KHwoGgswfP5elb7/Gsj1vKDD4RHdPL9v2HGLTzoNs3tnIWzsPsvdA61HtIqEAUycUM21iKdMmFHPaxBKmVhUTi/rtP1ERGa389n+jy4G709YtA65ggMDgnGtM/WyJ67GfAT4+TPX52oKJZ7P07ddYvmcNHzzzSq/LGZXqD7WzcUcDb+44yKa3G9i65xDdPUdOfRwOBZg2oYQZk0uZMamE6ZNKqR5XRDCocQYi4l++CQxmVg6UAPvSNtUA7x3iYS4FWpxzy9LWV5nZz4DxJOaMeIzEGIYR9fyLeVWzCQdCbK7fTmP7IUrzSrwuaURzzrF7fwvrttWzYVs9G3Y0sH+AroXJ44qYWV3GGVPKOH1yKVOqijUIUURyjm8CA9D3EIT0IeCdJGaTHIrrgB8OsH4LcJtzbrOZTQV+D5xJoqvjCGZ2A3ADQHV19RBP6w+xcIy542axct86lu9dy+XTL/K6pBHFOcfO2mbWbD7Aum0HWL+t/qhbGfOiIWZNKePMqeXMnFrOGdVlFOaFPapYRCRz/BQY+jp2o2nro8DRf7alMbMi4P3A36Rvc85dk/J+h5ndDfzAzG53zjWktb0fuB9gwYIFQ3ygrn+cP/EcVu5bx7I9bygwZEBNfSur36rjjc11rNtaT2PLkXm2rCjKWdMrmXNaObOnVVBdVax5DURkRPJNYHDONSRvf0wfrVcFbB3CIT4OPOecOzCEtlsBA6YBDYO0zSnzJ56NLX+IdbVv0t7dQV445nVJOaWto5s3Nh9g1Vv7Wb2pjn31Rw5QLC+OMnf6GObOqGTu9ArGVxboVkYRGRV8ExiSngXSZx1aAPxyCPtexwAzS5rZu4EC59xvUlZPTC53nkSNvlYaK+aMitPYVL+NN2o2cOHk87wuydf6uhlWbKxlxZv72bC9/ohHMhfkhTl7RiXzzhjD2TMqmTimUAFBREYlvwWGu4AlZnamc26jmV1NYqDivQBmdifwQWCBc66jbyczmwFMBZ4e4JiTgb80s2edc21mVkLicdyPOuf2D+tP45EFE89hU/02lu15Q4FhAD29cdZvref1DTW8tr7miDkQAgZnTi3nvFljOfeMMcyYVKq7F0RE8FlgcM6tMLNPkXh4Vd9Mj1emTNoUIzEAMv1PvOuAHzvn4hztOeA84A9m1gEUAc8DXx+GH8EXzp90Dj9d8ytW7ltHT7yXUEAPEOro6mHVpv38ce0+lq2vobWjp39bSWGE+bPGsWDWOObNHENRfsTDSkVE/MlXgQEg2XXwm2NsuwW4ZYD1tx3neLuAL2SswBwwoWgcE4uq2NNcw5t1mzlr3CyvS/JER2cPyzbU8uIbe1jx5n66ug/fRTt5XBEXzKnigjlVnF5dpoGKIiKD8F1gkMw4f9I57NlYw7I9a0ZVYOju6WX5xlpeXL2X1zfU0Nl1OCTMrC5j0dzxLJo7ngljCj2sUkQk9ygwjFALJpzNYxufZtmeN7ju3I+N6IF6zjk2bG/gDyt28dIbe2lt7+7fNnNKGRfPm8g75k5gTFmeh1WKiOQ2BYYRakbFVMrySjjQ1sBb9duYWTnd65Iyru5gO88u28nzy3dSU3944OK0CSW867yJXHTORMaWD3XOLxEROR4FhhEqYAEunrKQ37z5DC9sf3XEBIbunjivr6/h96+9zaq39uOSd0BWlMRYfN4kLp0/mSnji70tUkRkBFJgGMHeNfVCfvPmM/xx1wquO+/jRIK5O0Vx3cF2nn51B0+/9jaNzYnZFkPBAIvmjufyhdWcc/oYDVwUERlGCgwj2OSSCUwrq2bbwZ0s37OGd1TP97qkE+KcY83mAzzx8jZeX19DPHk1obqqiCsvnMLi8yZTXKBbIEVEskGBYYR719QL2XZwJy/seDVnAkN3Ty9LV+3hsRe2smNfEwChoHHR3Alc/c7TmH1a+YgexCki4kcKDCPcO6ecz4/feJQ3ajb4/pHXre3dPPnH7Tz+4jYOJrsdSouiXP2O07jqwimUFeu5GCIiXlFgGOGKo4WcN/4slu15gxffXsb7Z13udUlHaWzu5DcvbuW3L2+nLTkD49TxxXzgkum867yJhEOaqVJExGsKDKPAu6ZeyLI9b/DCjle5Zua7fXM5/2BTB7/4w2aeeuXt/lkYz55RyUcvO515Z4zxTZ0iIqLAMCqcN/4siiIF7Dy0hx2NuzmtbLKn9TS1dvHLP2zm8Ze29weFhbOr+Ni7T2fW1HJPaxMRkYEpMIwCoWCId045n6c2L+GF7a94FhjaOrr51ZKt/HrpVto7E10PF8yp4pNXzmLaRP+OrRAREQWGUWPx1At5avMSXtq5jE/P+0hWn2DZ2xvnmdd38tOn3+yfQ+G8WWP51JWzOKO6LGt1iIjIyVNgGCVOK6tmcskEdh3ay2u7V/LO6vOzct7lG2v578fXs6u2GUg82+HPr5nDnGkVWTm/iIhkhgLDKGFmvPf0xdy//CGe3PT8sAeGmvpWfvDYOl7fUAPAuPJ8Pvu+2Vx0zgQNZhQRyUEKDKPIxVMu4KE1v2Zzww7eOrCNMyqnZfwcXd29/HLJFn7+7Ft09cTJi4b4xBUzef/Fp+n2SBGRHKbAMIpEQxEun34Rj218mt++9XzGA8OaLXX8v5+/wb4DrQAsPm8Sf/7+OZRrwiURkZynwDDKXDVjMY+/+Qyv7V7FgdYGKgtO/TbGto5uHnhiA0+9sgNIPOvh/3z4bOZOrzzlY4uIiD8oMIwy5fmlLJo8n5d2LuOpLUv49DkfPqXjLd9Yy70/X82BQx2EgsYnrpjJRy47nVAwkKGKRUTEDxQYRqGrz7iMl3Yu47mtL/HR2VcTC594l0FHZw//9Zt1PP3q2wCcUV3KF/7kXKZUFWe6XBER8QEFhlFoRsVUZlZMY1P9NpbseJWrTl98Qvtv2d3Id36ygj11LYRDAT591Zl84F3TCQZ094OIyEil68aj1NUzLwPgd2/9gbiLD2mfeNzxyz9s5ivfW8qeuhaqq4r4l5vfxYcvnaGwICIywukKwyi1cOI8KvPL2deyn2V73uCCSecet31LWxfffWglyzfWAnDNO0/juvfPIRrWrZIiIqOBrjCMUsFAkGtnXQHAI2sfJx4/9lWGbXsO8aV7XmD5xlqK8sPc/pcXcOOHz1ZYEBEZRRQYRrF3T3snFfll7Graxx93rRiwzfPLd/GV7y2lpr6N6ZNK+NcvLWbh7KrsFioiIp7zXWAws2vNbJmZLTWzl81swSDt3zSzJWmvL6a1KTazB5PHXWlmd5vZqO+OCQfDfHT21QD8fP0T9MZ7+7f1xh3/9et1/Ov/rqSrJ87l51dz9+cvZlx5vlflioiIh3z1pWlm84GHgIXOuQ1mdg3wtJnNcc7VHGO3Gufc4kEO/SDQ6pw738wiwBLgG8DfZ6by3PWu0xbx2Ju/Z1/zfl58+3UWn7aIjs4evvPTFby2voZQ0LjxQ2dz5YVT9AwIEZFRzG9XGG4FnnbObQBwzj0B1AI3newBzews4EPAPyWP2QXcA9xsZoWnWnCuCwWCfGzO+wD4+frfUtfYwq33vcxr62soyAvzjRvewVWLpiosiIiMcn4LDJcDy9PWLQOuOMVjdgDr0o6ZB1x0CscdMS6qPp+JRVXUtdbzNz/+KVt2NVJVkc8///XFzJ2h6Z1FRMRHgcHMyoESYF/aphrgeE9JKjCz/06OeVhiZreZWerUhdOAWuecSztm37ZRLxAIcPH4SwHoKNvIzCnFfOcLlzB5XJHHlYmIiF/4JjAABcllZ9r6TuB4I+02Afc55y4BPgpcDfwi7bgDHZOBjmtmN5jZcjNbXldXN9Tac9r6bfU8/PMm4m1FWKSTCy/roKQw6nVZIiLiI34KDK3JZfo3VRRoO9ZOzrlPO+eWJd8fAG4H3mdm81KOO9AxGei4zrn7nXMLnHMLxowZc2I/QQ5atWk/d9z/Cu2dvcwMvxOAX7/5FAfaGjyuTERE/MQ3gcE51wA0Auk3+VcBW0/gUH1tZySX24CxduSovaq0tqPS6+tr+MYPX6Oru5crFlbzjU9eywWTzqWzt4sfr37U6/JERMRHfBMYkp4F0uddWJBcfxQzm2tm16etnphc7kwunyExwHFO2jHbgZdPqdoctmrTfr79o2X09Ma55qLT+PzH5hEMGJ+d91GiwQiv7lrJmpqNXpcpIiI+4bfAcBdwpZmdCWBmVwPjgXuTn+80s3UpgxorgL81s4rk9iiJWzNfA1YAOOfWA78CvpJsEwa+CNzjnGvJ1g/mJxu21/PNB1/vDws3fHAugeTDoyoLyvnQ7KsAeGDlI/T09nhZqoiI+ISvAoNzbgXwKeDHZrYU+BpwZcqkTTESAxX7uhfWkBjg+KSZLSFxxaAGuNY515ty6OsAzGwZiTDxR+COYf1hfGrr7ka+/l+v0tnVy7vPn8xffWDuUXMsvH/m5YwvHMue5hqe3Py8R5WKiIif2JF3G0qqBQsWuOXL06eFyF27apv56r0v0dTaxTvPnsBXPj2fYHDgzLh63wa+tfT7RENR7nnv/6UivyzL1YqIiAeOOUufr64wyPA52NTBHfe/QlNrF+fNGsuXP3XssAAwb/xsFk6cR2dPJ/cvfwgFSxGR0U2BYRTo7O7lzgde40BjO7OmlHHrZ88nHBr8H/1fnPcnFETyWbVvHc9teykLlYqIiF8pMIxw8bjjnv9dyVs7Gxlbns/X/vwCYpGhPXOsPL+U6+d/AoAfrX6UmpbRMZGViIgcTYFhhHvo92/y0ht7yYuGuOMvLqC06MRmcHxn9fm8o3oBnT2d3Pvqg8Tj8WGqVERE/EyBYQRbsmIXP3vmLQIGf/eZBUwZX3xSx7n+vE9QllfCpvpt/GbTMxmuUkREcoECwwi1fe8hvvfIagD+6oNzmT9r3EkfqzBawOfO/wwAP1v3ODsO7spEiSIikkMUGEagto5u7v7xMrp74lyxsJprLjr1h3LOGz+b98y4hN54L9/94w9o6WodfCcRERkxFBhGGOcc9/7iDfbUtTJ1fDE3fvjsjB37M+d8hNNKJ1PbUsf3X32QuNN4BhGR0UKBYYT5/Wtvs3TVHmKRIH/7ZwuIhoMZO3YkFOHLF91IYaSAVfvW8Yv1v83YsUVExN8UGEaQ7XsP8Z+/WgvATR89h8njijJ+jrEFFdy86C8xM36x/kmW71mT8XOIiIj/KDCMEB1dPf3jFt5zwRQWz588bOc6u+pM/nTuBwD4/msPsLe5dtjOJSIi/qDAMEL85HdvsqeuleqqIm740NxhP98HZr2HCyadS3t3B99eei+NHU3Dfk4REfGOAsMIsGF7Pb95cSuBgHHzJ87N6LiFYzEzblr4GaaVVVPbUse3l/4/2rs7hv28IiLiDQWGHNfZ3cv3frYK5+Ajl87g9MnZe6pkLBzj1ktuYlzhGLYf3MV3Xv5Penp7snZ+ERHJHgWGHPfQU4muiMnjivjT98zM+vlLYsV87V1/TUm0iLW1b3Lv6z/S7ZYiIiOQAkMO2/R2A4+9sIWAwc2fOJdwaPi7IgZSVTiGWy/5PLFQlJd3LueBlY/ocdgiIiOMAkOO6u7p5d9+toq4gw++awZnVGevK2Ig08qr+cpF/4dQIMTTW17ghyse1pUGEZERRIEhRz3+4jZ21bYwcUwBn7xqltflADB33Cz+9qL/QzgQ4vdbl3L/sp8qNIiIjBAKDDnoYFMHDz/zFgA3fPDsrNwVMVTzxs/h7y7+/4gEwzy//Y/8++s/1iOxRURGAAWGHPQ/v9tIe2cPC2dXcd6ssV6Xc5Szq87k7y/5PNFQlKU7XuOeV39IV2+312WJiMgpUGDIMZt3HeTZZTsJBY2/vHaO1+Uc0+yxZ/C1Sz5PXijGq7tW8o9L/o3mzhavyxIRkZOkwJBDnHP84LF1OAfXXjydCWMKvS7puGaNmcE33v1lKvLK2HRgK7c9+8/UNO/3uiwRETkJCgw55IVVe9i4o4HSwih/csUZXpczJFNKJ/HNy/+WqaWT2Neyn689+0+8WbfV67JEROQEKTDkiI7OHh58Yj0An7n6TPJjYY8rGrry/FK+ftmXOXf8WTR3tfL1Jf/K05tf0FwNIiI5RIEhRzz5x+3UH+pg+qQS3n1+tdflnLC8cIy/vej/cPXpl9Ib7+WHKx/m+68+QEdPp9eliYjIEPguMJjZtWa2zMyWmtnLZrbgOG3Hmdm3zOwlM1tiZqvM7FYzC6W1a0xuT339yfD/NJnR0dnDo3/YAsCfvfdMAgHzuKKTEwwEue68j/PFRX9BNBTlpZ3L+Ptn7mZPU43XpYmIyCBCgzfJHjObDzwELHTObTCza4CnzWyOc26gb5WPAVcCi51zzWY2GVgJ5AF3pLRb7ZxbPMzlD5sn/7idptYuZlaXcd5M/91GeaLeWX0+U0on8d2X72d30z6++vtv85l5H+Xy6RdhlpthSERkpPPbFYZbgaedcxsAnHNPALXATcdovx/4Z+dcc7L9LuAR4FNZqDUrUq8ufPLKWSPmC3VS8Xi+ffnfcdGUhXT2dvGDFQ9x14v3crD9kNeliYjIAPwWGC4HlqetWwZcMVBj59wjzrmH01a3A9FhqM0T/VcXppRx7swxXpeTUbFwjC9c+OfcvOgvKYjks2rfer781D/yx50rNCBSRMRnfBMYzKwcKAH2pW2qAaadwKEWkbjKkKrKzH6WHBfxrJl93sz8M5/yMbSnXl14z8i5upDuHdUL+O5VtzOvajYtXa3c88p/cfdL91HXWu91aSIikuSnMQwFyWX6sPlOIH8oBzCzy4Fq4H1pm7YAtznnNpvZVOD3wJkM0NVhZjcANwBUV3t7N8KTL4/cqwvpyvNKufWSz/Ps1pf4yZpfsnLvWtbXbuKjc97H+2a+m1DA9/lORGRE880VBqA1uUzvTogCbYPtbGZTgPuAa51zjanbnHPXOOc2J9/vAO4GPpe8qkFa2/udcwuccwvGjPHuS7q9s4dfLhl5YxeOx8y4YsbF3PPef+Adk+fT2dvFT9f8ir97+pusqdnodXkiIqOabwKDc64BaASq0jZVAcedGtDMxgG/Bq53zq0awum2AsaJdXVk1TOvv3346sIZI/vqQrqyvBJufsf1/P0lf824gkp2Ne3jzhe+x7eX/j92HdrrdXkiIqOSbwJD0rNA+rwLC5LrB2RmZcATwK3OuReS625I2f5uM7s2bbeJyeXOU654GMTjjide2g7ARy6dMSquLgxk3vjZfPeq2/nk2R8kLxRj1b713PL0ndy//CEa2hq9Lk9EZFTxW2C4C7jSzM4EMLOrgfHAvcnPd5rZOjOLJT8XAr8jESjqzGxBcqKnG1OOORn4ipnlJ/cpAb4EPOqc8+WTkFZu2s++A62MLctj4ez0Cy6jSyQU4YNnXsn33vd13jP9Egzj2a0v8te/vZ0HVj6i2zBFRLLET4Mecc6tMLNPAT82s3YgCFyZMmlTjMQAyL4/ub8IXJB8ffUYh30OOA/4g5l1AEXA88DXh+enOHWPv7gNgPe98zSCQb9lOm+UxIq5fsGfctUZi3lk7RO8unslv9v8B57d9hJXTL+Ya2deQXl+qddlioiMWKb73Y9twYIFbvny9Gkhhtfu/c187u7niYSDPHjHeyjKj2T1/Lni7cbd/Hzdb3l9z2ogMe30xdULef+sy5lcMsHb4kREctcx+8B9dYVB4LfJsQuXzp+ksHAcU0oncctFN7Lj4C5+ufEpXtu9iiU7XmHJjlc4d/xZXH3GpcwdN4uA6QqNiEgmKDD4SGt7N88tT4zDvOYi397A4StTyybzN+/4K2pa6nhi07Ms2f4Kq/atY9W+dYwvGst7pl/C4tMWURAZ0lQeIiJyDOqSOI5sd0n8ZulWfvDrdcydXsm3/r93Zu28I0lTZwvPbn2RZ7a8SH37QQCiwQiLJs/n0mmLmFU5eu86EREZAnVJ+F3qrZTvv/g0j6vJXcXRQj48+718YNZ7WLF3LU9veYG1tW/2d1dUFY5h8WmLuGTKBVQWHDVvl4iIHIMCg0+s3LSfffW6lTJTgoEgCyfNY+Gkeexr3s+S7a/wwo5XqWmp4+G1v+Hhtb9hVuV03ll9Posmn0dxrMjrkkVEfE1dEseRzS6Ju360jJfX7OUzV5/Jx959RlbOOdrE43HW1G5kyfZXWL53DV293QAELMBZY2eycNI8zp94DmV5JR5XKiLiGXVJ+FlLezevb6jBDC6dP9nrckasQCDAvPFzmDd+Du3dHSzfs4aXdi5jTc0G1tRuZE3tRn644mHOqDiN8yedw3kT5jKxqEpjHkREUGDwhZff2Et3T5yzZ1RSWZrndTmjQl44xsVTF3Lx1IU0d7awYu9aXtu9ijU1G9lUv41N9dv4yRu/YmxBBeeNn8u5E+Zw5pjTiYXSn40mIjI6KDD4wB9W7AIScy9I9hVFC1l82iIWn7aI9u4OVu1bz8q9a1lVs579rfU8tWUJT21ZQigQYmblNM4edyZnV53JaaWTCQQ0z4OIjA4aw3Ac2RjDUNvQxvXffIZIKMD/fP0q8mPhYT2fDF08HmdLww5W7lvLG/s2su3gThyH/3vJC8c4c8zpzBlzBrPHns7U0kkEA0EPKxYROWUaw+BXL6zcDcCFZ41XWPCZQCDAGZXTOKNyGp+Y+wGaO1tYt38Tb9RsZH3tJmpbD7By71pW7l0LQDQU5YyKqcysnMGsyulML5+iCaNEZMRQYPCQc+5wd8QCDXb0u6JoIYsmz2fR5PkAHGhtYEPdZtbvf4sNdZupbaljbe0m1tZu6t9nYlEVMyqmMqN8KtPLp1BdOpFIUMFQRHKPAoOHtuxuZPf+FkoKI5x7xhivy5ETVFlQziUFF3DJ1AsAaOxoYtOBrWyq28qm+m1sP7iLPc017Gmu4YUdrwIQtACTSyYwrayaqWWTmVo6ierSieSHNdhVRPxNgcFDf1iR6I645NxJeoz1CFAaK+aCSedywaRzAejp7WFH4262NOxgc/12th3cyd6mWnY07mZH427YfnjfcQWVVJdOZHLJBKpLJlJdMoGqorGENCZCRHxCgcEjPb1xlq5KBAbdHTEyhYKhRHdExVSuOn0xAB3dHexo3M22gzvZ0bibtxt3s+vQPmpbD1DbeoBle97o3z8YCFJVOIZJxeOZWFzFxKIqJhSPY3zRWF2REJGsU2DwyOq36jjU0sWksYXMmFTqdTmSJbFwjFljZjBrzIz+dT3xXvY21bDz0F52HtrDrkN72XloL3Wt9expqmFPU81RxymLlTC+aCxVhWOo6lsWjmVsYYXChIgMCwUGj7y4eg8Ai+dP0kyCo1woEKS6dCLVpROB8/vXd/R0sreplj1NNexu2sfe5lr2NtdS07yfgx2HONhxiA11m486XlG0kHEFlYwtqGBM8tX3viK/TJNPichJUWDwQG/csXxjLQCLzhrvcTXiV7FQlGnl1Uwrrz5ifTwe50BbAzUtdexr3p9Ytuxnf0uiW6O5s4Xmzha2NOwY8LhFkQIq88upKCinIq+UivwyKvLKqMgvpTwv8YqEIln4CUUklygweOCttw/S1NpFVUU+k8fpKYlyYgKBAGMLKxlbWMnZVWcesS3u4hzqaKa25QD7Ww9Q11qfeLXVs7+1gfq2gzR3tdLc1cr2xl3HPEdBJJ/yvFLKYiWU5hVTFiuhLK+E0lgJpbFiSvOKKY0VkxeK6QqZyCihwOCB1zck+qQXztaDjSSzAhagLC/x5T5rzPSjtsddnKaOZg60HeRAWyJA1Lc30tB2kPq2gzS0N9LQcYjWrjZau9rYdWjvcc8XDoYpiRZREiuiJFpEcayI4mjyfbSQ4lghRZFCiqOFFEULiYWi+ndeJEcpMHggNTCIZFPAApTmlVCaV8KMiqkDtom7OC2drTS0NybGSrQ3cTD5vrGjiUPtTTR2NtPYfojO3i4OtDVwoK1hSOcPBUIURQoojBZQHC2kMFJAYaSAomgBhZF8CiMFFETyKQjnJ7flUxDJV9AQ8QEFhiyrqW9lZ00z+bEQs6dVeF2OyFECFkhcKYgVMZXjz0Da0dPJoY4mDnU0c6izmaa+ZWcLTR3NNHe10NTZQnNnK02dzXT1dvcP2DwRZkZBOBEe8sMxCsL55IfzyI/kJZZHvGLkh/PI61uGYuSFY0RDEQKm+U5ETpYCQ5a9vj5xdeG8mWMJh/Q/L8ltsVCUWOEYxhUObabSzp4umrtaaOlMjKNo7myhpauNlq7W/mVfd0hLctna3UZHT2eyTetJ12pYot5wNBEiQrH+94n1yeUxXtFQJLmMEg1GEp+DUT2xVEYNBYYs6+uOuGCOuiNk9ImGIkRD5VTml5/Qfj3xXtq62mjtbqetu53WrraUZQdtyfVt3e20d3fQ3tPev76ju5O2ng46ezpp7+mgvaeDg5zYFY7jCQdCRJLhIRIKJ5cRosFIchkmEuz7nHgfCYaTr+O9DxFOfg4Hw0QCYYKBoLpmxDMKDFnU2t7Nuq31BALG/DPHeV2OSM4IBYL93SQnKx6P097TQUcyOHR0JwNEd2Jd4pV8392Zsq6Tzt5OOnq66OzpoqOnI7Hs7aKrp4vueA/dXT200pbBn3hgZkYkkAgQ4WCISCBMKLnsWxcOhAgFw4QDfe9DR78Phgn1rQsE+z8n3odStiVfwWD/9v6lJd4rxIwevgsMZnYtcDvQDgSBLzrnlh+nfTHwPWBOsv0zwNeccz0pbcYD/wGMBaLAQ8657wzbD3EMKzftpzfumDOtgqJ83ecukk2BQCAxoDKDjxx3ztHV201XbxedvYlA0dnTSVdvd8rnLrp6u/rXJdp305Wyvu8Yfe+7+9bFD3/u7u2m18UTx+3tytjPkAnBviBhgf4Q0Rcu+t9bkGAg2P85GAgdua7/fYCQBQn0tbOUfSyxvW9d0AL9+6SuC/S3S25P3ZZsGzBLLPva9b8/3E5jXo7kq8BgZvOBh4CFzrkNZnYN8LSZzXHOHT0/bsKDQKtz7nwziwBLgG8Af588ZgB4HPidc+52MysBVppZk3Pu/mH+kY6guyNERhYzS3azRMjGjCq98d5EeIj3JF693XT39tDV201PPLHsjifW9S371vfEe5LrE+t6+t/3JpbJtqnbe+K9/et6470DfO6h18XpjffSG++lMwu/g2xLhIcgAesLFokgkRo8ApbyPjDAur59LEAwEMD6Px/eFkxpk77NjjiGHdXusmnvpDRWPOy/C18FBuBW4Gnn3AYA59wTZlYL3ETiqsMRzOws4EPA2cn2XWZ2D/CgmX3LOdcCXA3MAy5PtjlkZv8J3GZmP3DOueH/saC3N86K5OyOC+eoO0JETlzfX+MxrwtJ4ZxLBIe+QOF6E2HD9R4OFf1ho5felPWHP8cT75OfU98fXhcnnjx+PB5PWx5e3+vixJPHiLs4vfFEoOlxic+JNonjxp1LtHG9/cEn8TnxHki87417/Fs+vvkT5o7KwHA5cHfaumXAFQwQGJLtO4B1ae3zgIuAp5JttjrnGtPaTAZmAm9movDBbNzRQHNbNxMqC5g0VrM7isjIYGbJ8RNhr0vJuHg8nggZfSHCJUNGMmz0h5C+V/xw276wcsT2I17uiM9Htj1ym8MNuL5vXUk0O98pvgkMZlYOlAD70jbVAO89xm7TgNq0qwQ1Kdv6lgMds29bVgLD6xv6ri6oO0JEJBcEAgECaBxDHz/9JgqSy/RusE7gWKOUCo7RnpR9htKmn5ndYGbLzWx5XV3doEUP1aGWTgKmwCAiIrnJN1cYgL4ZWdKfvRuFY96v1HqM9qTs0wqkd+6kt+mXHAh5P8CCBQsyNr7hS396Hn957VkUxPz0KxcRERka31xhcM41AI1A+p/gVcDWY+y2DRhrR94E3Lf/1pQ2Ax0ztU1WFBdECAZ98ysXEREZMr99ez0LLEhbtyC5fiDPkBjgOCetfTvwckqbGWZWmtZml3Nu06kWLCIiMhr4LTDcBVxpZmcCmNnVwHjg3uTnO81snZnFAJxz64FfAV9Jbg8DXwTuSd5SCfA7YDXwpWSbYuAG4M4s/UwiIiI5z1cd6s65FWb2KeDHZtY30+OVKZM2xUgMVEztgrgO+L6ZLUu2fxa4I+WY8eTskf9hZq8kj3F/tidtEhERyWWWpXmLctKCBQvc8uXHnJVaRERkpDnmg0H81iUhIiIiPqTAICIiIoNSYBAREZFBKTCIiIjIoBQYREREZFAKDCIiIjIoBQYREREZlAKDiIiIDEoTNx2HmdUBb2fwkJXAgQweb7TS7/HU6Xd46vQ7PHX6HZ66TP8ODzjnrhpogwJDFpnZcudc+sO15ATp93jq9Ds8dfodnjr9Dk9dNn+H6pIQERGRQSkwiIiIyKAUGLJLT8jMDP0eT51+h6dOv8NTp9/hqcva71BjGERERGRQusIgIiIig1JgyBIzu9bMlpnZUjN72cw0MvgEmVnEzO4ysx4zm+p1PbnEzK4xsyfN7Dkze9XMfmdmZ3tdVy4xs0vM7JdmtsTMXjSztWb2Ra/rymVm9nkzc2a22OtacomZXWdmbyb/XUx9lQzneUPDeXBJMLP5wEPAQufcBjO7BnjazOY452o8Li8nJAPC/wJvAUFvq8lJDwJfcM49BGBmdwHPmdlZzrlaTyvLHZ8EVjvnvgFgZucAK81sq3PuCW9Lyz1mNgH4itd15LC7nHMPZvOEusKQHbcCTzvnNgAk/+dSC9zkaVW5pRD4M+ABrwvJUUv7wkLSd0lM+PIej+rJRd8D/rXvg3PuDaARmOFVQTnu+8C3vC5Chk6BITsuB5anrVsGXOFBLTnJObfOObfF6zpylXPuw2mr2pPLaLZryVXOuQ3OuWYAMwuY2V8BncDPva0s95jZ+4Fu4Gmva5GhU2AYZmZWDpQA+9I21QDTsl+RCACLgA7gN14XkmvM7DZgL/Al4Grn3B6PS8opZlYAfJPE709O3jVm9ryZvWRmvzCz84f7hAoMw68guexMW98J5Ge5FhHMzIDbgducc/u9rifXOOfuBMaT+NJ7wcze4XFJueYfgf9wzqX/ESVDVwtsIRFYLwJ+BbxiZhcO50kVGIZfa3KZfuk3CrRluRYRSPQbv+2c+67XheQql/BTYClwl9f15AozOw+4APgPr2vJZc653znnvuqc60h+/inwCvDV4Tyv7pIYZs65BjNrBKrSNlUBW7NfkYxmZnYzMBv4iMel5BwzizjnutJWbwCu96KeHPU+IA94PnGhi1hy/T3J/0/e6Jzb5FFtuW4rMKy36+sKQ3Y8y9H/IBck14tkhZldD1wNfNw512Nm08zscq/ryiErkt05qSYAGsMwRM65f3TOneecW+ycWwx8Irnp5uQ6hYUhMLNvm1l6l/ZEYOdwnleBITvuAq40szMBzOxqEn2g93palYwaZvYJ4Gsk+t3nJicOuwK4yNPCcksR8Pm+D8n5VT4K/NCzimS0WgT8Zd8HM3sXcBnw78N5Uj1LIkvM7FoSA83aSUw8dLNzbpm3VeUOM4sAvwdKgXOA14C9A9wuKAMws24G7oL8unPuH7JcTk4ys0+S6H7IA3qTyx8C9zn9j/SEmdk9wIUkxjS8AWx2zn3M06JyhJldBfw1iRAbIPHf9necc78Y1vPq33MREREZjLokREREZFAKDCIiIjIoBQYREREZlAKDiIiIDEqBQURERAalwCAiIiKDUmAQERGRQSkwiIiIyKAUGETEd8zsm2b2Pa/rEJHDFBhEJGvM7J/N7OkhNH0f8ORw1yMiQ6fAICLZtBB4/XgNzGwCMANYko2CRGRo9CwJERl2yYeHtQDhlNUbnXOzB2h7PfBB59w1xzneJ4EvJY8ZA5Y4527NbNUikkqBQUSGnZkFgHOB5SSeTrgT6HTOHRyg7S+BZ5xz9x3jWBOS+5/hnNtmZmOAN51zFcP2A4iIuiREZPg55+LAeKAZWOacqzlGWIgAl3P88QvjSDwifmry2HXA1ZmuWUSOpMAgItlyLvCGO/5lzYuBXc65t4/TZjXwP8CzZva8md0ArMtcmSIyEAUGEcmWecCqQdpczSB3R7iEzwBzgRXAN4HVZlaagRpF5BgUGEQkW84B1gzS5mrgt8drYGYTzWyRc269c+4rwBxgAvDuzJQpIgNRYBCRbAkBs8xswkBXA8zsNBJf/C8NcpzTgbvNLJT8HAAM2JzBWkUkjQKDiGTL14BPALuBbw+w/X0k7o7oGeQ4bwJbgVfMbAnwOHCTc26wqxcicgp0W6WI+IKZPQn8wjn3317XIiJH0xUGEfGLJSSuFoiID+kKg4iIiAxKVxhERERkUAoMIiIiMigFBhERERmUAoOIiIgMSoFBREREBqXAICIiIoNSYBAREZFB/f+DXwfLHNWwHQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(time,y_analytical(time),label='$y$')\n",
    "plt.plot(time,dy_dt(y_analytical(time)),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 = 5.\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": "a4dd1a7ff1a14d2e8ddbfe6e56782669",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(BoundedFloatText(value=0.2, 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.2,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),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),ls=':',marker='d',label='absolute')\n",
    "    ax[0][1].plot(time_n[1:],(y_n[1:]-y_analytical(time_n[1:]))/y_analytical(time_n[1:]),\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$ and $\\Delta t = 0.3$?\n",
    "- Can you recover the forward and Backward Euler solutions?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}