04b_Generalized_Midpoint_Newton_Nonlinear.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 4 - Generalized midpoint rule, Newton 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": 13,
      "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",
        "%run plot_functions/plot_settings.py"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {},
      "outputs": [],
      "source": [
        "y_eq = 1.\n",
        "k = 2."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "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": 16,
      "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": 17,
      "metadata": {},
      "outputs": [],
      "source": [
        "time = np.linspace(0,5,100)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtQAAAIWCAYAAAB3DRpPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB6PklEQVR4nO3deXhU1f0/8PedNZkkM9kzZCFshk02MShhVcCVoq1YBQSB7qCItGoR+hVXqlWrCAp1CagUqaK4/KiWoiJKZQlJkCXsCSH7OpNl9rm/PyYZMiSBLJPcmcz79Tz3mTvn3rn3M5LSdw7nniOIoiiCiIiIiIg6RCZ1AURERERE/oyBmoiIiIioExioiYiIiIg6gYGaiIiIiKgTGKiJiIiIiDqBgZqIiIiIqBMYqImIiIiIOoGBmoiIiIioExioJSCKIoxGI7imDhEREZH/Y6CWQE1NDXQ6HWpqaqQuhYiIiIg6iYGaiIiIiKgTGKiJiIiIiDqBgZqIiIiIqBMYqImIiIiIOoGBmoiIiIioExioiYiIiIg6gYGaiIiIiKgTFFIXQERERERtZ7PZ4HA4pC7DL8nlciiVSq9fl4GaiIiIyA8YjUaUl5fDYrFIXYpfU6vViI6Ohlar9do1GaiJiIiIfJzRaERBQQFCQ0MRHR0NpVIJQRCkLsuviKIIm80Gg8GAgoICAPBaqGagJiIiIvJx5eXlCA0NRWJiIoN0JwQHByMsLAwXLlxAeXm51wI1H0okIiIi8mE2mw0WiwU6nY5h2gsEQYBOp4PFYoHNZvPKNRmoiYiIiHxY4wOIXfEwXaBq/G/prYc7/SpQWywWbN26FdOnT4der0dUVBRiYmJw++2347///W+7r/fee+8hNTUVsbGx0Ov1uPvuu3Hy5MlWz8/JycHMmTOh1+sRGxuL1NRUbN68uTNfiYiIiKhN2DvtPd7+b+lXgfqRRx7Bvffei+HDh+PUqVOoqKhARkYGrFYrpk2bhrVr17b5WitXrsT8+fOxePFilJSUICcnB3a7HampqTh8+HCz8zMzMzFmzBg4HA7k5OSgpKQEixcvxrx587Bq1SovfksiIiIi8ieCKIqi1EW01QMPPIDDhw/ju+++82gvLy9H79694XQ6UVxcjPDw8MteJyMjA6mpqZg7dy42bdrkbjcYDEhKSkJKSgoOHDjg/u1FFEWMGjUKZ8+eRX5+PnQ6nfsz8+bNw+bNm5GRkYGRI0e26XsYjUbodDoYDAavTtlCREREPY/ZbMa5c+fQt29fBAUFSV1Oj+Dt/6Z+1UN9yy23YMWKFc3ao6OjMWjQIFgsFhw6dOiK11m3bh1EUcSsWbM82nU6HW699VZkZGRg79697vY9e/YgOzsbt912m0eYBoBZs2bB6XRi3bp1HfxWREREROTP/CpQT58+HTfffHOLx6xWKwAgKirqitfZtWsXAGD06NHNjjW27dy5s8PnExEREVHg6BHzUJeXl+PUqVMYMmQIhg8fftlz6+vrcf78eahUKsTExDQ7npCQAMD1AGKj48ePexxrKiYmBkqlEnl5eTCZTAgODu7MV+kSHx3dgXxDISx2C/48cbHU5RARERH1KD0iUK9ZswZ2ux1r1qy54lOb1dXVAACNRtPi8ZCQEABAVVVVs880HmtKEARoNBoYDAZUV1e3GKgtFovHMqFGo/GyNXpbdvExnCg/AwCwOWxQyjntDhEREZG3+NWQj5b8+OOP+Otf/4qnnnoKU6ZMkbqcFq1evRo6nc69JSUldev9w9Sh7v0aa1233puIiIjIG5599lkIgtDqlMWbNm2CIAhYvXp1N1fm54H62LFjmD59OpYsWdLiw4otaZwBpL6+vsXjdXWuwBkREdHsM43HmhJF0X2t1mYXWb58OQwGg3vLz89vU63eolVd7FmvtTBQExERkf8ZO3YsAGD//v3NjtXU1ODPf/4z+vfvj2XLlnV3af475OPIkSOYOnUqFi5ciBdeeKHNn9NoNOjduzfOnz+PsrKyZuOoCwoKAACDBg1ytw0ePNjjWFNlZWWw2WxITk5udfy0Wq2GWq1uc43e1rSH2miplawOIiIi8r6H//4tqmosVz5RQhFhavz94cmdusaYMWMgl8uxb9++ZseefPJJFBcX49NPP5Ukc/lloM7MzMRNN92ExYsXeyyqkpubC5VKhfj4+Mt+fsqUKUhPT0dGRgZuueUWj2MZGRkAgGnTpnmc/9RTT7mPXel8X+M55IOBmoiIqCepqrGgwmCWuowuFxoaimHDhiErKwtWqxUqlQoAcOLECaxZswY33XQTZsyYIUltfheoDxw4gJtvvhl//vOf8eijj3ocW7VqFfr06eMO2aIo4sKFC83GLC9evBgbN27Eli1bPAK1wWDAv//9b1xzzTVIS0tzt0+YMAEjRozAjh07YDAYPOai3rJlC2QyGRYv9t3ZM8KaDPmoYQ81ERFRjxIRJt2/greVt2pMS0tDVlYWsrOzkZqaCgBYunQpRFHEK6+84pV7dIRfBeq9e/fi1ltvRXx8POrr65st+Z2VlYU+ffq43y9ZsgRr167FkiVL8Oqrr7rbR48ejccffxyrV6/GjTfeiHnz5sFoNGL+/PkAgHfeecdjthBBEJCeno6JEydiwYIFSE9Ph1arxaZNm7B582asXLmyzaskSsFzyAfHUBMREfUknR1K4U/S0tLw+uuvY9++fUhNTcVnn32GL7/8EkuXLnUP0ZWCXwXqF154AUajEUajEU8++WSL59x5553u/aSkJGg0mhZn1XjmmWeQkpKCV199FY888ggEQcCECRNw4MABDBw4sNn5o0aNwoEDB7BixQqkpKTA6XQiOTkZmzZtwn333ee179gVtE0CdS17qImIiMhPNX0w0WKxYNmyZYiJicETTzwhaV2CKIqipBUEIKPRCJ1OB4PBAK1W2+X3K6opxUM7XD9o45PHYMn1C7r8nkREROQdZrMZ586dQ9++fREUFCR1OZLT6/XQ6XS4//77sWLFCvzjH//Ab37zm3Zdw9v/Tf162jxqmzA1x1ATERFRzzB27FicOnUKzz77LEaNGoVf/epXUpfEQB0INMpg95hwBmoiIiLyZ2lpae51QNasWQOZTPo4K30F1OVkgsw90wcDNREREfmzvn37AgBmzZqF8ePHS1yNCwN1gGic6cPIpceJiIjIT4miiL///e/QaDTtWtivqzFQB4jGmT4sdgusDpvE1RARERG13xtvvIG9e/fiiSeeQGJiotTluPnVtHnUcWGqplPn1SFSEy5dMURERERtdOHCBbz//vs4efIk3nvvPYwfPx5//OMfpS7LAwN1gPBc3KWWgZqIiIj8wn/+8x8sX74ckZGRuPfee/Hqq69CLpdLXZYHBuoA4TF1npUPJhIREZF/WLhwIRYuXCh1GZfFMdQBoumQD870QUREROQ9DNQBwnNxF870QUREROQtDNQBQttkDDWHfBARERF5DwN1gLj0oUQiIiIi8g4G6gDRNFBzDDURERGR9zBQB4jGpccBjqEmIiIi8iYG6gChUQZDJrj+uDmGmoiIiMh7GKgDhCAI7mEf7KEmIiIi8h4G6gCibRj2wTHURERERN7DQB1AGnuoLQ4rLHarxNUQERER9QwM1AEklMuPExEREXkdA3UA0XosP85x1EREROTfzp8/j5ycHKnLYKAOJJyLmoiIiHqSefPmYfDgwVKXwUAdSMK4/DgRERGR1ymkLoC6j1bNIR9ERETUc3z77bdSlwCAPdQBJazJQ4lGDvkgIiIi8goG6gASpuIYaiIiIiJvY6AOIGEe0+ZxyAcRERGRNzBQBxDO8kFERET+6tlnn4UgCNi8eXOLxzdt2gRBELB69epuroyBOqAEK4Igl8kBMFATERGRfxk7diwAYP/+/c2O1dTU4M9//jP69++PZcuWdXdpnOUjkAiCAK0qFFVmA2f5ICIiIr8yZswYyOVy7Nu3r9mxJ598EsXFxfj000+hVqu7vTb2UAeYxuXHjZyHmoiIiPxIaGgohg0bhqysLFitVnf7iRMnsGbNGtx0002YMWOGJLWxhzrANM5FbXPYYLFboVaoJK6IiIiIOuPP/1mNarNR6jIuKzxIi7/etLzT10lLS0NWVhays7ORmpoKAFi6dClEUcQrr7zS6et3FAN1gLl06jy1IlLCaoiIiKizqs1GVJqqpS6jW6SlpeH111/Hvn37kJqais8++wxffvklli5dKukS5AzUAebSxV2iQxioiYiI/Fl4kFbqEq7IWzU2fTDRYrFg2bJliImJwRNPPOGV63cUA3WA8Zg6j+OoiYiI/J43hlL4i379+iEuLg779u3DSy+9hDNnzuAf//gHwsPDJa2LDyUGGC3noiYiIiI/NnbsWJw6dQrPPvssRo0ahV/96ldSl8RAHWhCVZ5DPoiIiIj8SVpaGkRRRH19PdasWQOZTPo4K30F1K2a9lDXcvlxIiIi8jN9+/YFAMyaNQvjx4+XuBoXvw3UR44cQVpaGgRBQG5ubrs+O3nyZOh0Ouj1+mZbXFwcBEHA4sWLm30mMjKyxc/ccccdXvxmXavpGGr2UBMREZE/EUURf//736HRaPDCCy9IXY6b3z2UaDab8fTTT2P9+vWd6uJ/9dVXMX/+/GbtO3bswO23345Zs2Y1O/bxxx9j8uTJHb6nL/B4KJGrJRIREZEfeeONN7B37148//zzSExMlLocN78L1EuXLkVxcTGys7Nx3333Yffu3e2+xpgxY5CQkNDisfXr12P48OE+808I3qZtMoaaDyUSERGRr7tw4QLef/99nDx5Eu+99x7Gjx+PP/7xj1KX5cHvAvXy5cuRnJzcqWu09k8E+fn52LFjB9auXdup6/sytUINhUwBu9POQE1EREQ+7z//+Q+WL1+OyMhI3HvvvXj11Vchl8ulLsuD3wXqzobpy3nzzTcREhKC++67r8vuITVBEBCmDkGVyYAaPpRIREREPm7hwoVYuHCh1GVclt8+lOhtdrsdb7/9NubOnYvQ0NAWz9m6dSvGjx+PpKQk6PV6TJo0CRs2bIDD4ejmajtH27D8eI2lFqIoSlwNERERkX9joG7w+eefo7CwEIsWLWr1nNzcXGzYsAH5+fnIzs7GDTfcgEWLFmH69Omw2Wytfs5iscBoNHpsUmp8MNHmtMNit0haCxEREZG/Y6BusH79ekyaNAlDhgxp8fiHH36IL774AkOHDgUAxMXFYdWqVZg9eza+/PJLvPHGG61ee/Xq1dDpdO4tKSmpS75DW3lMncdhH0RERESdwkAN4OzZs9i5c+dle6djYmJaHAB/9913AwC2bdvW6meXL18Og8Hg3vLz8ztfdCeEqTnTBxEREZG3+N1DiV1hw4YN0Ov1+PnPf97uz8bHxwMAioqKWj1HrVZDrVZ3uD5vC1NxLmoiIiIibwn4Hmqr1Yr09HT8+te/hlKpbPGcrKwsvPnmmy0eKywsBADExsZ2WY3exh5qIiIi/8OJBLzH2/8te3SgFkXxisMrtm3bhsrKSvzud79r9ZysrCw89thjMJvNzY598sknAIAZM2Z0rthupG26WqKVgZqIiMiXNQ45vdwECNQ+jf8tvTWfdY8O1EuWLEHv3r3x0EMPtXrO+vXrMWPGjFZXTmxUVVWFuXPnoqCgAABQW1uL559/Hu+++y7GjRuHJUuWeLX2ruTxUCJ7qImIiHyaUqmEWq2GwWBgL7UXiKIIg8EAtVrd6uiE9vK7MdT79+939wZXVlYCAFJTUyGXyzFnzhy89NJL7nOTkpKg0WhanVUjJycH3333HXbu3HnZe951110QBAHbtm3DxIkTUVNTA7PZjIEDB+KFF17Agw8+CJVK5aVv2PU8x1AzUBMREfm66OhoFBQU4MKFC9DpdFAqlRAEQeqy/IooirDZbDAYDKitrb1iZ2p7CCJ/1el2RqMROp0OBoMBWq222+9fVleBxV+sBABcn3gNlo37TbfXQERERO1jNBpRXl4Oi4VrSHSGWq1GdHS0VzOY3/VQU+eFcQw1ERGR39FqtdBqtbDZbH63SrOvkMvlXhvm0RQDdQBSy1VQypWwOWycNo+IiMjPKJXKLgmF1HE9+qFEapkgCNA2jKPmGGoiIiKizmGgDlCNc1EbrbV8YpiIiIioExioA1TjOGqH0wGTvfn82kRERETUNgzUASpMdXG1xFqOoyYiIiLqMAbqAMXFXYiIiIi8g4E6QHHqPCIiIiLvYKAOUNqmgZpDPoiIiIg6jIE6QDXO8gFwyAcRERFRZzBQB6gwVdMeagZqIiIioo5ioA5QnmOoOeSDiIiIqKMYqANU0yEf7KEmIiIi6jgG6gClVYe59w1mo4SVEBEREfk3BuoApZIrEaLSAACqGKiJiIiIOoyBOoBFBukAAJWmaoiiKHE1RERERP6JgTqARQSHAwBsDhvqbPXSFkNERETkpxioA1hEsM69X2UySFgJERERkf9ioA5gDNREREREncdAHcAiG4Z8AK5x1ERERETUfgzUAaxpoGYPNREREVHHMFAHMA75ICIiIuo8BuoA1jRQV5qrpSuEiIiIyI8xUAew8KAmPdT11dIVQkREROTHGKgDmEImh65hCfJKM4d8EBEREXUEA3WAaxz2UW0ywCk6Ja6GiIiIyP8wUAe4xtUSHaITNZZaaYshIiIi8kMM1AGOU+cRERERdQ4DdYDzmOmDgZqIiIio3RioA1ykx1zU1dIVQkREROSnGKgDXASXHyciIiLqFAbqABcRxNUSiYiIiDqDgTrARXqslshATURERNReDNQBTqsOg0xw/RhwDDURERFR+zFQBziZTOYe9sEhH0RERETtx0BN7qnzDOYaOJwOiashIiIi8i8M1OQO1CJEGMw1EldDRERE5F8YqOmSxV2qpSuEiIiIyA/5baA+cuQI0tLSIAgCcnNz2/XZ3NxcyOVy6PX6Frd//vOfLX4uJycHM2fOhF6vR2xsLFJTU7F582YvfBtpRXIuaiIiIqIOU0hdQHuZzWY8/fTTWL9+PWSyjv8+kJSU1K4gnpmZiUmTJmHKlCnIycmBTqfDpk2bMG/ePJw6dQqrVq3qcC1S41zURERERB3ndz3US5cuxdGjR5GdnY2hQ4d2yz1FUcSCBQsAABs3bkR4eDgEQcD8+fMxZ84cPP3008jKyuqWWrpC09USq8zVktVBRERE5I/8LlAvX74c27dvR2JiYrfdc8+ePcjOzsZtt90GnU7ncWzWrFlwOp1Yt25dt9XjbR6Lu7CHmoiIiKhd/C5QJycnd/s9d+3aBQAYPXp0s2ONbTt37uzWmryp6RhqDvkgIiIiah+/C9TeUl9fj2XLlmHo0KGIi4tDv379MHv2bBw4cKDZucePHwcAJCQkNDsWExMDpVKJvLw8mEymLq+7K4SoNFDKXMPpGaiJiIiI2idgA3VVVRX0ej2+//57FBYWYvv27cjPz8fYsWOxceNGj3Orq6sBACEhIc2uIwgCNBqNx3mXslgsMBqNHpsvEQTBPXUelx8nIiIiap+ADNRJSUkoKirCo48+ioiICMjlcgwfPhyffvopQkNDsWjRIpSUlHjtfqtXr4ZOp3NvSUlJXru2tzQ+mFhjrYPVYZO2GCIiIiI/EpCBWi6XIzo6ull7ZGQkbrzxRphMJuzYscPdHh4eDgCoq6tr9hlRFFFfX+9x3qWWL18Og8Hg3vLz8zv/Jbys6eIu1Rz2QURERNRmfjcPdVeLj48HABQVFbnbBg8eDAAoKChodn5ZWRlsNhuSk5MRHBzc4jXVajXUanUXVOs9kUGeM33Ehjb/hYOIiIiImgvIHuqNGze2+PAhABQWFgIAYmNj3W1TpkwBAGRkZDQ7v7Ft2rRp3i6zW0Vqwt37nIuaiIiIqO16dKAWRbHF4RUbN25Eenp6s/bq6mp8++23UKlUuOWWW9ztEyZMwIgRI7Bjxw4YDJ7DIbZs2QKZTIbFixd7/wt0o4igcPc+Z/ogIiIiarseHaiXLFmC3r1746GHHmp27M0338Rbb70Fq9UKADhz5gxmzpyJ6upq/O1vf/NYOEYQBKSnp7tXTDQYDBBFERs3bsTmzZuxcuVKjBw5sru+VpeI4OIuRERERB3id4F6//790Ov10Ov12Lt3LwAgNTUVer0ef/zjHz3OTUpKgkajaTarxoYNG7By5Uq89dZb6Nu3LyIjI5GWlobQ0FDs2rULS5YsaXbfUaNG4cCBAxAEASkpKYiNjcXatWuxadMmPPnkk133hbtJ09USOXUeERERUdsJoiiKUhcRaIxGI3Q6HQwGA7RardTlAABMNjPu//hhAMDQ2BQ8ccPDEldERERE5B/8roeaukawMghBCtdMJBxDTURERNR2DNTkdnG1RAZqIiIiorZioCa3yIbVEk12M0w2s7TFEBEREfkJBmpya1x+HACqzOylJiIiImoLBmpy85zpg4GaiIiIqC0YqMktIohT5xERERG1FwM1uTUd8lHJQE1ERETUJgzU5BbJ1RKJiIiI2o2BmtwiOIaaiIiIqN0YqMmNY6iJiIiI2o+BmtxUChVCVSEA2ENNRERE1FYM1OShcdhHpdkAURQlroaIiIjI9zFQk4fGBxNtDhvqbPUSV0NERETk+xioyUNEULh7n8M+iIiIiK6MgZo8RHhMnVctXSFEREREfoKBmjx4BOr6aukKISIiIvITDNTkIVoT4d4vr6+UsBIiIiIi/8BATR5iQ6Ld+6V1FRJWQkREROQfGKjJQ2xIlHu/tK5cwkqIiIiI/AMDNXkIUgZBpw4DAJTUMlATERERXQkDNTXT2EtdZTLA6rBJXA0RERGRb2OgpmZiQ13jqEWIKOc4aiIiIqLLYqCmZpo+mFjCcdREREREl8VATc3EhTaZ6aOWPdREREREl8NATc2wh5qIiIio7RioqZnYpj3UDNREREREl8VATc1EBYdDLrh+NEo5dR4RERHRZTFQUzNymRzRDVPnldSVQxRFiSsiIiIi8l0M1NSixrmoTTYz6qz1EldDRERE5LsYqKlFcXwwkYiIiKhNGKipRU0fTOQS5EREREStY6CmFjWdOo8zfRARERG1joGaWuS5uAsDNREREVFrGKipRRxDTURERNQ2DNTUohCVBsHKIABAaR2XHyciIiJqDQM1tUgQBHcvdXldBRxOh8QVEREREfkmBmpqVeODiQ7RiQpTtbTFEBEREfkoBmpqVSwfTCQiIiK6Ir8N1EeOHEFaWhoEQUBubm67PnvhwgU888wzuPbaaxEVFYXw8HAMGDAADzzwAIqKilr8zOTJkxEZGQm9Xt9su+OOO7zwjXxPHKfOIyIiIroihdQFtJfZbMbTTz+N9evXQybr2O8DQ4YMQVBQEN5//31MnToVAPDVV1/hvvvuw4cffoi9e/eif//+zT738ccfY/LkyZ0p36/Ehka597m4CxEREVHL/K6HeunSpTh69Ciys7MxdOjQDl3D6XRi9erVuOmmmyCTySCTyXDrrbfiqaeeQmlpKVasWOHlqv0Te6iJiIiIrszveqiXL1+O5OTkTl3jsccew89+9rNm7RMmTAAA7N27t1PX7ymiQ6IgQIAIkWOoiYiIiFrhd4G6s2EaAP7yl7+02G61WgEAUVFRLR4PNCq5EhHBOlSaqjkXNREREVEr/G7IR1dq7JmeNWtWi8e3bt2K8ePHIykpCXq9HpMmTcKGDRvgcFx+jmaLxQKj0eix+YvGJcgNlhqYbWaJqyEiIiLyPQzUDcxmM9atW4dBgwbhgQceaPGc3NxcbNiwAfn5+cjOzsYNN9yARYsWYfr06bDZbK1ee/Xq1dDpdO4tKSmpq76G18V6jKNmLzURERHRpRioGzz66KMoKyvDhx9+CI1G0+z4hx9+iC+++ML9IGRcXBxWrVqF2bNn48svv8Qbb7zR6rWXL18Og8Hg3vLz87vse3hbbMjF4S98MJGIiIioOQZqAM8//zzS09OxY8cOXH311S2eExMTA7lc3qz97rvvBgBs27at1eur1WpotVqPzV/Ehca49zl1HhEREVFzAR+o//a3v2H16tX46quvcP3117f78/Hx8QDQ6oIw/q7pkI8S9lATERERNRPQgfqZZ57B888/j127diEtLc3dfvDgQfeMHwCQlZWFN998s8VrFBYWAgBiY2O7tliJNF3chWOoiYiIiJrr0YFaFMVWxyuvXLkSr732Gr755huMHj3a41hqaqo7KAOuQP3YY4/BbG4+y8Unn3wCAJgxY4YXK/cd4UFaKOVKAOBc1EREREQt8Lt5qNtjyZIlWLt2LZYsWYJXX33V3f7II4/gxRdfxL333ott27Zddvxzo6qqKsydOxevvPIKEhISUFtbi3Xr1uHdd9/FuHHjsGTJkq78KpKRCTLEhkShwFiM0rpyiKIIQRCkLouIiIjIZ/hdoN6/f7+7N7iyshKAq0dZLpdjzpw5eOmll9znJiUlQaPReExTV11djRdffBEA8MEHH7TpnnfddRcEQcC2bdswceJE1NTUwGw2Y+DAgXjhhRfw4IMPQqVSeesr+pzYkGgUGIthddhgMBsRHqyTuiQiIiIinyGIoihKXUSgMRqN0Ol0MBgMfjHjxzsZW/Hl6W8BAE9P+RMGRveXtiAiIiIiH9Kjx1CTdzR9MJFT5xERERF5YqCmK+JqiUREREStY6CmK4oLbRKo2UNNRERE5MGrgbqqqgrvvvuuNy9JPsCzh5qBmoiIiKgprwbq8+fPY8GCBd68JPmAYGUQwtShALhaIhEREdGl2jVt3vnz5y97vOliKNSzxIVEo8ZSi8r6atgcNvdiL0RERESBrl2Buk+fPpdd1IOLfvRcsaHROF2ZCxEiSurKkajtJXVJRERERD6hXYE6IiICTz/9NCZNmtTi8ZycHPzyl7/0SmHkWxK1evf+BUMRAzURERFRg3YF6tGjR6OqqgpDhw5t8bjdbgfXiemZmgboC8ZiCSshIiIi8i3tCtR/+MMfUFdX1+rx3r17Iz09vdNFke9J1DUJ1AaOlSciIiJq1K5A/fOf//yyxyMiInD//fd3qiDyTfrQWMhlcjicDvZQExERETXBhV2oTRQyOeJDYwEAhTUlcDgdEldERERE5Bs6HajlcjlOnjzpjVrIxyXq4gEAdqcdJbVlEldDRERE5Bs6Haj5EGLgaDrTR76xSMJKiIiIiHwHh3xQm3k+mMhATURERAQwUFM7JGnj3fvsoSYiIiJyYaCmNtOHxkAuuH5kCthDTURERASAgZraQSFXoFdYHACggDN9EBEREQFgoKZ2alwx0e60o6SuXOJqiIiIiKTHQE3twgcTiYiIiDx1OlA/9thjiIqK8kYt5Acae6gB4AIfTCQiIiJq39LjLVm9erU36iA/kdSkhzrfUChhJURERES+gUM+qF16hca6Z/q4YCyWuBoiIiIi6TFQU7so5Arow2IBAIXGYjidTokrIiIiIpKW1wL1qVOn8NRTT+H11193t33xxRfo168fwsPDsXTpUi5T3kM0jqO2caYPIiIiIu8F6nfeeQcvv/wyKioqAAAlJSWYNWsWCgsL0b9/f6xfvx5vvvmmt25HEuKDiUREREQXeS1Q79y5Ex999BH+8pe/AADeffdd1NXV4cMPP0RGRga2bt2Kt99+21u3IwnxwUQiIiKii7wWqEtKSjB16lT3++3btyMlJQU/+9nPAAB33HEHLly44K3bkYQ8e6j5YCIREREFNq8F6qbjo0tKSrBv3z7MnDnT4xxBELx1O5JQr7BYyBpn+mAPNREREQU4rwXq2NhYfP/99wCAN954A6Io4q677nIfP3/+PEJDQ711O5KQUq6EPjQGAFBQU8KZPoiIiCigdXphl0YLFy7ErbfeioEDByIzMxPXXXcdRo0aBQDYv38/Hn/8cYwZM8ZbtyOJJep6obCmBDaHDaV15e6p9IiIiIgCjdd6qBcvXoyHHnoIJpMJt9xyCzZv3gwAcDqduOeee3DmzBnMnj3bW7cjiXGmDyIiIiIXQeTk0N3OaDRCp9PBYDBAq9VKXU6H/HD+AF793zsAgFnD7sDPh9wicUVERERE0uBKidQh7KEmIiIicmGgpg6JD4tzz9pywcBATURERIGLgZo6pOlMHxdqijnTBxEREQUsBmrqsCRtPAC4Zvqor5C4GiIiIiJp+G2gPnLkCNLS0iAIAnJzczt0jffeew+pqamIjY2FXq/H3XffjZMnT7Z6fk5ODmbOnAm9Xo/Y2Fikpqa6ZzMJRIk6vXufwz6IiIgoUHU4UK9fvx6lpaXerKVNzGYzVqxYgUmTJuHUqVMdvs7KlSsxf/58LF68GCUlJcjJyYHdbkdqaioOHz7c7PzMzEyMGTMGDocDOTk5KCkpweLFizFv3jysWrWqE9/If/HBRCIiIqJOBOpFixYhISEBkyZNwtq1a1FQUODNulq1dOlSHD16FNnZ2Rg6dGiHrpGRkYHnnnsO9913H+bPnw9BEBAeHo6NGzdCFEUsXLjQYyl1URSxYMECAMDGjRsRHh4OQRAwf/58zJkzB08//TSysrK88fX8SmLDkA+APdREREQUuDocqPPy8vDCCy/A4XDgoYceQnJyMsaOHYuXX365w0Mw2mL58uXYvn07EhMTO3yNdevWQRRFzJo1y6Ndp9Ph1ltvRUZGBvbu3etu37NnD7Kzs3HbbbdBp9N5fGbWrFlwOp1Yt25dh+vxV/HaizN9nDd0zy9URERERL6mw4E6KSkJDz/8ML7//ntcuHABr776KoKDg/HYY4+hf//+GD16NP76179edkxyRyQnJ3f6Grt27QIAjB49utmxxradO3d2+PxAoZIrkRjmGkedbyiE1WGTuCIiIiKi7ueVhxJ79eqFxYsX4+uvv0ZRURHWr1+PmJgYPPHEExg8eDCGDRuGp556CkeOHPHG7Tqlvr4e58+fh0qlQkxMTLPjCQkJAFwPIDY6fvy4x7GmYmJioFQqkZeXB5PJ1EVV+65+ka5fcByiE3nVFySuhoiIiKj7eX2Wj+joaPzmN7/Bl19+iZKSErzzzjvo27cvVq9ejREjRmDQoEH44osvvH3bNquurgYAaDSaFo+HhIQAAKqqqpp9pvFYU4IguK/VeN6lLBYLjEajx9ZT9I+8+C8GZyrzJKyEiIiISBpdOm1eeHg47r//fnz22WcoKyvD+++/j6uvvhonTpzoytv6nNWrV0On07m3pKQkqUvymgGRfdz7DNREREQUiBTddaPQ0FDMmjWr2YOA3S08PByAa+hHS+rq6gAAERERzT7TeKwpURTd12o871LLly/HsmXL3O+NRmOPCdXJ4QmQy+RwOB0M1ERERBSQ/HZhl47SaDTo3bs3rFYrysrKmh1vnP5v0KBB7rbBgwd7HGuqrKwMNpsNycnJCA4ObvGearUaWq3WY+splHIleutc0+cVGIthspklroiIiIioewVcoAaAKVOmAHDNR32pxrZp06Z1+PxA079h2IcIEeeqzktbDBEREVE369GBWhRF5OfnN2tfvHgxBEHAli1bPNoNBgP+/e9/45prrkFaWpq7fcKECRgxYgR27NgBg8Hg8ZktW7ZAJpNh8eLFXfMl/MAAjwcTGaiJiIgosPToQL1kyRL07t0bDz30kEf76NGj8fjjj+P999/Hpk2bIIoiDAYD5s+fDwB455133AuWAK6ZPNLT090rJhoMBoiiiI0bN2Lz5s1YuXIlRo4c2Y3fzLd4zvSRK10hRERERBLwu0C9f/9+6PV66PV692qGqamp0Ov1+OMf/+hxblJSEjQaTYsPAD7zzDNIT0/HmjVrEBcXh5SUFMjlchw4cAAjRoxodv6oUaNw4MABCIKAlJQUxMbGYu3atdi0aROefPLJrvmyfiJR2wsquRIAZ/ogIiKiwCOIoih29iLPP/88li1bBqVS6Y2aejyj0QidTgeDwdBjHlD8y3//hhMVZwEAb9/5N4SpQyWuiIiIiKh7eKWHevny5bj66qvx//7f//PG5cgPNR32cZYPJhIREVEA8UqgTk9PR01NDWbMmIHbbrsNJ0+e9MZlyY/05wIvREREFKC8Eqjvv/9+nDx5Eg8//DB27dqFYcOG4U9/+lOPWmKbLq9/ZG/3/mkGaiIiIgogbQ7UJ06cgEwmg0LR8uKKoaGhePHFF3H48GHceOONePnll5GSkoJ33nnHa8WS79KHxSJYGQSAM30QERFRYGlzoNZoNJg4cSImTJhw2fMGDhyIf//739i+fTtCQ0Px61//Gqmpqfjf//7X6WLJd8kEGfpHuMZRV5kMqDRVS1sQERERUTdpc6BOSkrCt99+i2+++aZN58+YMQPHjh3Ds88+i+PHj2P8+PGYO3cuCgsLO1ws+TaPBxM57IOIiIgChNfnobbb7cjIyMDrr7+O3/72t3j33XdhMpkgiiI2b96MgQMH4qWXXoIXZusjH9M0UHMcNREREQWKlgdEt9OWLVuwf/9+7Nu3D1lZWbBYLO7AHBcXhzvuuANpaWno27cvXn75ZTzyyCP4/PPPsX37doSHh3ujBPIBnismMlATERFRYPBKoJ4zZw4AQC6XY9iwYUhLS0NaWhrGjh2Lvn37epx711134fXXX8fSpUuxdOlSbNy40RslkA+I1kRCqw6F0VKLs5V5EEXRYwl3IiIiop7IK4H6mWeewdixYzFmzBiEhIRc8fxFixbh4MGD+Pzzz71xe/IRgiCgf2QfZBYdQY21DmV1FYgNjZa6LCIiIpKYw+GE1e6E1eaA1eaEze5o8t61b2t8tbvOsdodsNmczdr6J4Tj5uuTr3zTbuSVQP3444+3+zMpKSmorKz0xu3Jh/SP7I3MoiMAXOOoGaiJiIh8hyiKsDtEd5C1NAZam7PJfpP39ovn2TzOcbqPNdu3XQzLNrsDFpsTTqf3np0bPyK+Zwbqjpg7dy7i4uKkuj11Ec8VE3OR1nu0dMUQERH5CadThNXugMXqCq+Xf7W7Aq/1YiBuad/SLDS7jnkx20rCZndKXUIzkgXqhIQELFiwQKrbUxfxmDqv6ryElRAREXmP0ynCbLW7w6rZ6nC9tzS8Nm2zuoKr2X3uxbam781WzyDcE6gUMiiVcqiVMqiUcigVrn2lQg5VQ5tKIYdSKYNKIYdK0XBe4/uGc5WKxnMbryNzf06rUUn9NZuRLFBTzxQepEWUJgIV9VU4W3keTtEJmeD12RmJiIha5HCKsFjtMFnsMFsdrtcW9s3WhlfLJa9We/NjVv8MvHKZALVK7gqmSjnUTYLuxfeNx2XutovHPc9tDLvqhgDsPrdJ6A3UyQgYqMnr+kcko6K+Cia7GYU1JUjU9pK6JCIi8lE2uxMmi/3iZm54tTbZt7hCrslsR32TfXdIttphtjhQb7H7RfCVyQQEqeQIUsmhVioaQq8MQSqFK7yqGsJui68KqJUy977KvS/3+KxKKYdCzg6t7sJATV7XPzIZ+wuyAABnKvIYqImIehiHU4TJYke92QaT2Y46sw315ouB12RpeG+xo97ccF7DvsnScE7Dvt3he+NhATQEXgWC1K5XdUMAvrivcAXiS/bVjfvKi593h2GVwtW7q2DQ7WkYqMnrBkT1ce+fqDiLSX2vl64YIiLyYLM7UGdyheA6kw31ZhvqzHbUmxpeG8Jx42vdpaHZYoPJ4ju9wCqFDEFqBYLUCmjUrrAbrFYgWO0KtkFqBYJVF/dd7xv23aFZ4RGg1Up5wA5doI5hoCavuyqqL+QyORxOB46WnpC6HCKiHsVqc6DOZEOtyebx2hiQPdrc7Q2B2GSDVeIZEpQKmTvwaoIU7n331tCmaQzFlx5XKxCkcp0XpOKwBvIN3RaoDx8+jAEDBkCj0XTXLUkiQQo1BkT2wYnyMyiqKUWlqRqRweFSl0VE5DNsdgdq6m2orbeipt4VfGvqrag12VDb0F5rdu27ArLVvS9VIG4MwK5NCY1aAU1ww2uQ0n0sWK10B+WLr0p3GOZwB+qJui1Qjxw5EsHBwViwYAFefPFFBAUFddetSQJDY6/CifIzAIBjpScxPnmMxBUREXmfze5ETb0VNXVWGOutqK23wljXGJRdYbmmYb+2/mJotli7d8iEJkiBkGAlQoKUCAl2Bd5L9zXBSoQEXQzHIUFKBDe+qhWQyTgEgqg13RaoZ82ahfPnz2PDhg04dOgQdu/eDaVS2V23p242NHYgPj72JQDgaOkpBmoi8nk2uwPGOqtrq214rbO4Xutd7xuDc03DeeZuCsbBajlCglUIDXaF4NZeQ5q81zQE5mC1AnKGYaIu1W2BevPmzQCA06dPY+bMmVi/fj0efPDB7ro9dbOUqH4cR01EkjJb7KiudQViQ62lYbPC0PC+MTAbGsKzyWLv0noUchm0IUqEalQI07jCcahGeXE/WIkQjQphGmXDsYsBmuOEiXxbtz+UOGDAAOzcuRNarba7b03dSK1QISWqL46XnUZxbRkq6qsQpYmQuiwi8mMOp4iaOiuqay2orjGjutaK6hrXvqHW1W6otbgDc1cNq5AJQFiIKxSHaVTQNu6HuMJws/2GYKxWceYIop7KK4H6+eefx7Jly9o8hCMmJsYbtyUfNyQmBcfLTgMAjpaexMQ+10lcERH5GlEUUWuyocpoRlWNBVUNAbnKaEFVjbnhvcXV01xrgVP0fg1hGiW0IWpoQ1TQharc+41bWJN9rUYFTZCS44mJyIMgimKn/3qSyWS46qqr8PLLL+P222/3Rl09mtFohE6ng8Fg6NE99UdKTuCpb18BANzYNw2/HzNX2oKIqNs4nSIMdRZUGS2oNJpRYTCjqsaMSqPZFZ6NFlQ2BGdvLuwhCGgIxmqEh6rd+7oQFbSNbaGucKwLUSNMo4ScwymIqJO80kOdnp6O5cuXY8aMGbj55pvxyiuvICUlxRuXJj+WEtUXCpkCdqcdR0tPSl0OEXmJyWJHhcGEimozKowmVBjMqDSYUWFseDWYUFVjgcNL3ckKuYDwUDXCw9QIDwtqsq+GLlSNiFA1dGGusBwWouIDeETU7bzSQw0AtbW1WLVqFV577TUAwIMPPoj/+7//69E9sB0VKD3UALDq65dxrOwUAOD16c8iOiRS4oqI6HLqTDaUV5tQbjChvNoVjsurGzaDGZUGE+rM3nl4TxuiQqQ2CBFhakRc+hoWhPAwNSLC1AgJVnLsMRH5NK8F6kYnTpzA0qVL8dVXXyE2NhbPPfccFi5c6M1b+L1ACtT/OvIFPjr6/wAAi8fcz2XIiSRkdzhRYTCjtKoeZVUmlFW7XsurTShrCM31XgjL4aFqRGqDEKkLcr1qgxCpdbVFNLwPD1Nz5goi6jG8HqgbffbZZ1i2bBnOnj2L0aNHY82aNRg7dmxX3MrvBFKgPlZ6Equ++TsAYHLfsVg0Zp7EFRH1XFabA6VV9SitNLleL9mvNJrRmb/xVUo5onRBrk0bjOhwV2iO0gUjqiFAR4QFcSU8Igo4XTZt3owZM3DLLbfgpZdewrPPPovx48dj9uzZeP755xEfH99VtyUfMyCqL5QyBWwcR03UaQ6niIpqE4or61BcUY+SynqUVNSjpLKuITBbOnxtpUKG6PBgxIQHI7px0wVd3A8PRiiHXhARtcjrgdputyM7Oxv79u3D/v37sW/fPphMJoiiiM2bN2P79u1YtWoVli1bxr+YA4BKrkRKdD8cLT2JsroKlNZVIDYkSuqyiHyW2WJHcWU9isrrUFReh+KKhq2yHmVV9bA7OtbFHB6qRkxEsGsL1yC2yX50eDB0oSr+nUxE1EFeCdRbtmxxh+esrCxYLBY0jiSJi4vDHXfcgbS0NPTt2xcvv/wyHnnkEXz++efYvn07wsPDvVEC+bChsSnu3uljpScR25dDfyiwmS12FFXUobC8DoVltSgqd+0Xldd2uJc5IkyN2EgN4iI0iIvSICbCtR8bGYyYCA3USrmXvwURETXySqCeM2cOAEAul2PYsGFIS0tDWloaxo4di759+3qce9ddd+H111/H0qVLsXTpUmzcuNEbJZAPGxp7cQrFo6UnMZmBmgKAw+FESWU9CspqUVDmCs4FDVuFwdzu6wWr5YiLDIE+SgN9VAj0kRrERYUgLlKDmIhgBKm6feFbIiJq4JW/gZ955hmMHTsWY8aMQUhIyBXPX7RoEQ4ePIjPP//cG7cnHzcgsg+UciVsDhuOlp6EKIr8p2XqMUwWOy6U1iC/pBYXSmtwobQWF0prUVRe2+7hGeGhavSKdoXmXlEhrv3oEPSKCoE2hEMyiIh8lVcC9eOPP97uz6SkpKCystIbtycfp5QrMSi6H34qOYHy+kqU1VUgNjRa6rKI2qXOZMP54hqcLzEiv6QW+SU1yC+tQVmVqV3X0YaoEB8dgviYUMTHhCA+KhS9YkIQHx0CTZCyi6onIqKuJNm/Ec6dOxdxcXFS3Z662ZCYFPxUcgIAcKT0JG5koCYfZbE5kF9cg9wiI/KKjThfUoPzRUaUt2OYhkIuQ0JMCBJiQ5EYG4aEGFeATogJRZhG1YXVExGRFCQL1AkJCViwYIFUt6duNjR2IADXEJ9jpSdxY780aQuigCeKIkoq63Gu0OgKz0VG5BYZUFReh7aumB0SpEBSXBiS4sKQGBuGpDhXgI6N1HD5ayKiAOKXT7Hk5ORg5cqV+P777+F0OpGcnIylS5e6H468ksmTJyMzMxPBwcHNjomiiNLSUixatAjr1q3z+Mzhw4ehUjXvXbruuuvw6aefdvwLBYABkclQy1WwOKwcR03dzmZ3IK+oBmcLDThXYMC5IiNyCw1tXkI7JFiJ3nFhSO6lRVJcKJLjtEjShyEiTM2fYyIi8r9AnZmZiUmTJmHKlCnIycmBTqfDpk2bMG/ePJw6dQqrVq1q03VeffVVzJ8/v1n7jh07cPvtt2PWrFnNjn388ceYPHly575AgFLIFRgY3R+HS46jwlSF4toy9AqLlbos6oHMFjvOFhpwtsCAMxcMOFNQjfPFNXC0odtZpZAhSR+GZL0WfXppkdxLi2R9GCK1QQzORETUKr8K1KIouoeJbNy4ETqdDgAwf/58fP3113j66adx5513YuTIkZe9zpgxY5CQkNDisfXr12P48OEYP368V2snYFjcIBwuOQ4AyCj8CdMHTpG4IvJ3VpsD5woNOJ1fjVMXqnE6vxr5JTVtGrIRHR6MvvFa9IvXoU+8K0D3ig7lUA0iImo3vwrUe/bsQXZ2Nu655x53mG40a9YsvPfee1i3bh3efPPNy17nhRdeaLE9Pz8fO3bswNq1a71WM100OmEYNh/+BACQUXiYgZraxekUUVBWixN5VTiZX4WT56uQW2i8Ys+zTAAS48LQP0GHfg1b33gdHw4kIiKv8atAvWvXLgDA6NGjmx1rbNu5c2eHr//mm28iJCQE9913X4evQa1LCNNDHxqD4toyHC87jVpLHULVV563nAJTbb0VOXlVyMmrxIlcV4iuv8KYZ7lMQLJei/6JOgxICkf/BB2Se2m56AkREXUpv/p/mePHXcMFWhquERMTA6VSiby8PJhMphYfOLwcu92Ot99+G3PnzkVoaGiL52zduhUrV65EXl4ebDYbBg4ciNmzZ+PXv/415HIu63slgiDg2vjh+OLkLjhFJzKLjmJCnzFSl0U+QBRFFJXX4di5Chw7V4mcvErkl9Re9jOCACTFhWFAYjhSksIxICkcfeJ1XGKbiIi6nV8F6urqagBocTVGQRCg0WhgMBhQXV3d7kD9+eefo7CwEIsWLWr1nNzcXGzYsAFDhw5FSUkJ3njjDSxatAjbt2/HZ599BqWy5UUZLBYLLBaL+73RaGxXbT3JtQmuQA0ABwsPM1AHKIfDibOFBhw7V4mjZytw/Fwlqmstl/1MpFaNlN4R7u2qpHAuhEJERD7BrwJ1V1q/fj0mTZqEIUOGtHj8ww8/RGRkpLsnOi4uDqtWrcKZM2fw/vvv44033sCSJUta/Ozq1avx5JNPdlnt/mRgdH+EqDSos9Yjq+go7A47FHL+GPZ0docTZy5U48iZCvx0phzHzlXCZGl9+IZcJqBfgg6D+kRiUHIEBiVHIiYimDNtEBGRT/KrJBMeHg4AqKura3ZMFEXU19d7nNdWZ8+exc6dO/HBBx+0ek5MTEyL7XfffTfef/99bNu2rdVAvXz5cixbtsz93mg0IikpqV019hRymRzX9Loae/L2w2Q341jZKQzXD5a6LPIyp1PE2UIDDp8qR/bpMhw/VwGTxdHq+SFBCgzqE4mh/aIwuE8kruodwaEbRETkN/wqUA8e7ApeBQUFzY6VlZXBZrMhOTm53cM9NmzYAL1ej5///Oftrik+Ph4AUFRU1Oo5arUaarW63dfuqa5NGI49efsBAAcLDjNQ9xCF5bXIOlmG7FNl+Ol0OWrqba2eGx6mxrD+0RjaNxJD+kWht17L6eqIiMhv+VWgnjJlCp566ilkZGQ0O9bYNm3atHZd02q1Ij09Hb///e9bHQOdlZWFAwcO4De/+U2zY4WFhQCA2FguUtJWI/RDIJfJ4XA6cLDwMBZc80v+U74fqjPZcPh0GTJPlCHzZCmKK+pbPTdSq8bV/aNxdf9oDOsfhYSYUP6ZExFRj+FXgXrChAkYMWIEduzYAYPB4DEX9ZYtWyCTybB48WJ3myiKuHDhwmWHV2zbtg2VlZX43e9+1+o5WVlZeOyxxzB37lwEBQV5HPvkE9e8yjNmzOjo1wo4GmUwro5NQXbxcZTXV+K8oQDJ4YlSl0VXIIoicouMOHCsBBk5JcjJq4KzlTmgQ4KVGD4gGiMGRGP4VTFIjGWAJiKinsuvArUgCEhPT8fEiROxYMECpKenQ6vVYtOmTdi8eTNWrlzpsUrikiVLsHbtWixZsgSvvvpqi9dcv349ZsyY0erKiY2qqqowd+5cvPLKK0hISEBtbS3WrVuHd999F+PGjWt1/DS1bHT8cGQXu6ZBPFhwmIHaR5ksdmSfKsPB4yU4eLwEFQZzi+cp5AKG9I3CyJQYjEyJQb+EcA7hICKigOFXgRoARo0ahQMHDmDFihVISUmB0+lEcnIyNm3a1GxBlqSkJGg0mlZ7qHNycvDdd99dcTGYu+66C4IgYNu2bZg4cSJqampgNpsxcOBAvPDCC3jwwQehUnHVtfa4Nn443jm0FYBr+ry7ht4mcUXUqNJoxv6jxdh3tBhZJ8tgdzhbPC8hJgSjBsbimoGxuLp/NILVfvfXCRERkVcIoiheft1e8jqj0QidTgeDwQCtVit1OZJ59KtnkVt9AQCwfsZqRAaHS1tQAMsvqcGPR4qw70gxTpyvavEclUKGYQOikTo4DqMHx0EfxVUuiYiIAD/soaae49qE4e5AfajwJ0ztP0HiigKHKIrIK67B3sOF+OFwIc4X17R4XpQuCGOG6pE6OA7DBkRzCW8iIqIW8P8dSTLXxg/HR0d3AHCNo2ag7lqNDxXuySrA3sOFKChrPp87APTppcV1Q/W47mo9BiSG82FCIiKiK2CgJsn0jeiNiGAdqkwG/FSSA7PdgiAF5+v2tsLyWuzJLMDuzALkl7TcEz24TyTShvfC9Vf34lAOIiKidmKgJskIgoBr44dj55k9sDntOFx8HGMSR0pdVo9QVWPGd5kF2H3oAk7lVzc7LgjA0H5RGDc8HmOH9UKUrn2LIREREdFFDNQkqWsTXIEacM32wUDdcVabA/uPFWPXgXwcOlHa4hzRg/tEYtKoBKQNj0eENqiFqxAREVF7MVCTpIbGDoRaoYbFbsHBgsOwO+xQyPlj2VaiKOJUfjX+u/88vssqQJ2p+XLf/eJ1mDgqARNGJiA2UiNBlURERD0bkwtJSiVXYnT8MOw9fxC11jpkFh9FasIIqcvyebX1Vnx76AK++jEPuUXGZsejdUG44dok3DA6CUlxYRJUSEREFDgYqElyE5Ovw97zBwEA3+XuY6BuhSiKOHauEl/9mIsfsgthtXsuuKJWyTF2WC9MuTYJwwbEcKVCIiKibsJATZIboR8MnToMBksNMgp/Qq21DqEqzjTRyGSx49uMfHzxw7kW54semByBm69LxrgR8dAEKSWokIiIKLAxUJPk5DI5xiWnYsfJr2F32vFj/iHOSQ2goKwWO344h/8eOI96s93jWGiwEjdcm4Sbr0tGcq/AXW2TiIjIFzBQk0+YmHwddpz8GoBr2EegBmpRFJF1sgzbvzuDQzmlzY4P7hOJ29L6YOzweKiVcgkqJCIioksxUJNP6BuRhERtL1wwFiGn/AxKa8sRGxotdVndxmZ3Yk9WAT759nSzhwxVChkmXZOI28f1Rf/EcGkKJCIiolYxUJNPEAQBE/tch38e3g4A+C5vP2YOvU3aorpBncmGr37MxWd7zqLCYPY4FhsRjNvH9cXUMcnQhqgkqpCIiIiuhIGafMb45FRsOfwpRIjYk7sPdw25FYLQM2eqMNRa8Ol3Z/D/fjjXbHx0Su9w/GLyVbh+WC/O1EFEROQHGKjJZ0RrIjEk9iocLT2JotpSnK7MxVVRfaUuy6uqjGZ8/O1p/Pt/ubBYHe52QQDGDNHj55MHYEjfyB77iwQREVFPxEBNPmVi8nU4WnoSALA798ceE6grDCZ89PUp/OfHPI/5oxVyATde2xs/n9wfibFcgIWIiMgfMVCTT7kuaRTeOvQBbA4b9p7PwPyRd/v1UuSGWgs++voU/t8P52BrEqSVChluvi4Zv7jhKsREBEtYIREREXWW/yYV6pE0ymCkJozw+6XI6802bN99Btt3n4bJcnFoh1olx61j++DnkwcgUhskYYVERETkLQzU5HP8eSlyq82BL74/h4++Pomaepu7XaWQYfr4fvjFDQOgC1VLWCERERF5GwM1+Rx/XIpcFEV8n1WIjf/vKEqrTO52uUzATdcn456pKYjScWgHERFRT8RATT7n0qXI/3f+EKYN8N2VE0+er8Jbnx7B8dxKd5sgAJOuScTsmwahV7Rv/zJAREREncNATT5pUp/r3UuR/+fMd5jaf7zPTSVXVmXCuzuO4dtDFzzaR6XEYOGMq9Gnl1aiyoiIiKg7MVCTT+obkYSrIvvgVGUu8qov4HjZaQyJvUrqsgC4lgnfvvs0Pth5ElbbxQcOE2ND8asZV2P0oFifC/9ERETUdRioyWfdmnIDTv2YDgD496lvfCJQ/3SmHG9sO4z8khp3W5hGidk3D8ItY/tAIZdJWB0RERFJgYGafNb1idfg3aBtqDYbsb8gC+V1lYgOiZSkluoaC9K/OIqvD+a722QCcFtaX8y5ZRBCNSpJ6iIiIiLpMVCTz1LIFbhpwET868gXEEURX53ejTkjft6tNYiiiP/sO4+NXxxFreniNHgpvcPxh7tGYEBieLfWQ0RERL6HgZp82tT+E/DxsS9hd9qx6+wPmDn0dqgV3dMbXFpVj7X/ykLmyTJ3W0iwEvffNhg3Xd8HchnHSRMREREDNfm48CAt0nqPxne5+1BrrcOevP2Y2n98l97T1Sudh7c/OwqTxe5unzw6EQt/NhQRYVzhkIiIiC5ioCafd9tVN+C73H0AXA8nTuk3rstm0WipVzpKF4QH7h6JawfHdck9iYiIyL8xUJPP6xeZjIFR/XCi4izyDYU4WnoSV8cN9Oo9RFHENxn5WP/xTx690tPG9MbCGVcjNFjp1fsRERFRz8FATX7h1pQbcOJ/ZwG4eqm9GajrzTa8se2wxwIt7JUmIiKitmKgJr8wJnEUIoJ1qDIZcLDwMEpryxEbGt3p6548X4UX389AUUWdu+3Ga5PwmzuHsVeaiIiI2oSrUJBfUMjkuHnAJACu4Rlfnt7dqes5nSI+/uYUHn1tjztMa4IU+NOc0Xh41jUM00RERNRmDNTkN6b2Gw+lzPWPKl+f/QFmm7lD1zHUWvDkWz8i/YtjcDhFAMDA3hF4ddlkTLom0Wv1EhERUWBgoCa/oQ0Kw7jkVABAvc2E/5zZ0+5rnL5QjYdf2Y1DJ0oBAIIA3D3lKvz1gfHQR4V4tV4iIiIKDAzU5Fd+NnAqBLimzPss5z8w2y1t/uzXB/Px2Gt7UFZlAgCEh6nx9G/TMO+2IVDI+T8FIiIi6himCPIrSbp4jE26BgBgtNTiq1NXHkttdzjx5vaf8Pcth2C1OwEAA5Mj8MrDkzAiJaZL6yUiIqKej4Ga/M7Mobdf7KU+sfOyY6mrayz4y4a9+GzPWXfbzdcnY/WicYjSBXd5rURERNTz+WWgzsnJwcyZM6HX6xEbG4vU1FRs3ry5zZ/Pzc2FXC6HXq9vcfvnP//ZJfcl70jU9UJa79EAgBpLbaszfpwvNmLZq7tx5EwFAEAhF7B45gg8cPdIKBXybquXiIiIeja/C9SZmZkYM2YMHA4HcnJyUFJSgsWLF2PevHlYtWpVm6+TlJSE4uLiFrfZs2d32X3JO2YOvd29/PhnOTthuqSX+qfT5Xi0yXjpSK0aqxeNxy1j+3R3qURERNTDCaIoilIX0VaiKGLUqFE4e/Ys8vPzodPp3MfmzZuHzZs3IyMjAyNHjrzsdXJzczF58mTk5uZ2630bGY1G6HQ6GAwGaLXaNn2GmlvzYzq+z9sPALh32Az8YsitAIDdhy7glQ8yYXe4xkv3T9Th/351PSK1QZLVSkRERD2XX/VQ79mzB9nZ2bjttts8Qi0AzJo1C06nE+vWresx96XLmzn0Nncv9ecn/os6az0++voUXtyc4Q7T1w6Ow+pF4xmmiYiIqMv41dLju3btAgCMHj262bHGtp07d/aY+9LlxYfFYULyGHyXuw911no8+ekHOP5jlPv4zdcn4w+/GA45p8QjIiKiLuRXSeP48eMAgISEhGbHYmJioFQqkZeXB5PJdMVr1dfXY9myZRg6dCji4uLQr18/zJ49GwcOHOjS+5J3zRxyG2SC68f4nDUTkNsAAHNvHYzFM0cwTBMREVGX86u0UV1dDQAICWm+op0gCNBoNB7nXU5VVRX0ej2+//57FBYWYvv27cjPz8fYsWOxceNGr97XYrHAaDR6bOQdkUFR0Fr7AgAEhR2qXnlYNvsa/HJqins4CBEREVFX8qtA7S1JSUkoKirCo48+ioiICMjlcgwfPhyffvopQkNDsWjRIpSUlHjtfqtXr4ZOp3NvSUlJXrt2IDNb7Xj6nX0oPpYAUXSFZ01iPkYN1V3hk0RERETe41eBOjw8HABQV1fX7Jgoiqivr/c4rzVyuRzR0dHN2iMjI3HjjTfCZDJhx44dXrvv8uXLYTAY3Ft+fv5l66MrM1nseOqtfcg6WQbRogEqEgEAFqcFWw5/KnF1REREFEj8KlAPHjwYAFBQUNDsWFlZGWw2G5KTkxEc3PEV8OLj4wEARUVFXruvWq2GVqv12Kjj6kw2PPGP/+GnM+UAAE2QAo9Ouw/BStdMHt+c24vTFbkSVkhERESBxK8C9ZQpUwAAGRkZzY41tk2bNu2K19m4cWOLDx8CQGFhIQAgNjbW6/elzqs12fCXDXtxPLcSABASrMTTv0tDakpv/HLodPd57xzaCqfolKpMIiIiCiB+FagnTJiAESNGYMeOHTAYDB7HtmzZAplMhsWLF7vbRFFscXjFxo0bkZ6e3qy9uroa3377LVQqFW655ZYO35e6htlqx1Nv/YhT+dUAgDCNCs/9YRxSekcAAG6+ajIStb0AAKcrc/Fd7j6pSiUiIqIA4leBWhAEpKenQxRFLFiwAAaDAaIoYuPGjdi8eTNWrlzpsVrhkiVL0Lt3bzz00EPNrvXmm2/irbfegtVqBQCcOXMGM2fORHV1Nf72t78hMTGxw/cl77PZnfjrpgPunmldqAqrF41Dv4SLDyAqZHLMH3W3+/3mw9tRb+VUhkRERNS1/CpQA8CoUaNw4MABCIKAlJQUxMbGYu3atdi0aROefPJJj3OTkpKg0WiazaqxYcMGrFy5Em+99Rb69u2LyMhIpKWlITQ0FLt27cKSJUs6dV/yLqdTxCsfHEJGTikA15jpJ38zFsm9mo9FH64fjDGJIwEABrMRHx3b0ewcIiIiIm8SRFEUpS4i0BiNRuh0OhgMBj6geAWiKGL9x4exY28uAEClkOHJ347F1f2bz9LSqLSuAg//+0nYHDbIBRlevOUvSNDqu6liIiIiCjR+10NNgWXzVznuMC2TCXjs/tTLhmkAiA2Jwh2DXA+JOkQnNmb+C/y9kYiIiLoKAzX5rM++O4OtO0+63z987yiMGdK2nuY7Bt2MKI3rYcXs4uM4UJDdJTUSERERMVCTT/rfT4V489Mj7ve/vXMYJo9u+wqTaoUK80be5X7/dsYHqLU0X5iHiIiIqLMYqMnnnC0w4KV/HnK/v3faQPxsQr92X+f6xGswUj8EAFBlNuCdzH95rUYiIiKiRgzU5FOqasx4+p19sFgdAIDJ1yRi9s0DO3QtQRDwu9T7oFG6VrD8Pm8/9l/I8lapRERERAAYqMmH2OwOrN54AOXVrrmjU3qH48FfjoQgCB2+ZpQmAgtG/dL9/h8HN8Norul0rURERESNGKjJJ4iiiHUfZbsXbonSBWHFguugUso7fe2Jfa7DtQkjAABGSy3ezNjCWT+IiIjIaxioySds330Guw64lolXKeVYueA6RGqDvHJtQRDw22tnI0wVAgDYdyETP5w/6JVrExERETFQk+QOHi9B+hdH3e8fnjUKA5LCvXqP8CAtfjV6lvv924c+QJXJ4NV7EBERUWBioCZJlVbW48XNGWgcgTHrpoEYPyKhS+6V1ns00pJGAwDqrPXYcOB9Dv0gIiKiTmOgJsnYHU787f2DqDPZAABjh/XCvdM6NqNHW/1q9L3QBbmWez9UdARfnd7dpfcjIiKino+BmiTzz69ykJNXBQCIjdRgyT2jIJN1fEaPtghTh+J3185xv9+U9RFOlp/t0nsSERFRz8ZATZLIPFGKj74+BQCQywQ8et9ohAYru+Xe1yYMx/SBUwEADqcDL+99EwazsVvuTURERD0PAzV1u6oaM17ecsg9bnrebYMxMDmyW2uYM/xODI65CgBQaarGq/97Bw6no1trICIiop6BgZq6ldMp4u//PITqGgsA4JpBsbhz0oBur0Muk2Pp2F8hvGE89ZHSE9h65PNur4OIiIj8HwM1dauPvz2NzJNlAICIMDUevveaLh833ZqIYB0eTvs1ZILrfwbbj3+FAwXZktRCRERE/ouBmrrNibxKvPfv4wAAQQD+OHs0wsPUktY0OOYq3DfiF+73a/dtRHFNqYQVERERkb9hoKZuYbU58MoHmXA6XQOnfzklBSNSYiSuyuX2lBtxfdI1AACTzYy//bAB9VaTxFURERGRv2Cgpm6x9b8ncaG0FgBwVVI4Zt3UtfNNt4cgCPhD6lwkhOkBAPmGQry0dwPsDrvElREREZE/YKCmLne2wIBtDVPkKeQCHrpnFORy3/rRC1YG4ZEJv0eYKgQA8FPJCbx+4D2upEhERERX5Fuphnoch8OJNf/KhKPJUI/kXlqJq2pZfFgcHpuwCEq5az7s7/P2Y8tPn0pcFREREfk6BmrqUtt3n8GZCwYAQLI+DDOnpEhc0eWlRPfDQ9cvhCC4Zh7ZfvwrfHWKy5MTERFR6xioqcsUlNXin1/lAABkArDknlFQKnz/R25M4kgsHHWP+/07h7Zi/4Us6QoiIiIin+b76Yb8ktMp4rV/ZcFqdwIAZkzsj5TeERJX1XY3XzUJdw6+GQAgQsSrP76DE+VnJK6KiIiIfBEDNXWJr37MxdGzFQAAfZQGc24eJHFF7Tdr2B2YkDwGAGBz2PDc7rU4WX5W4qqIiIjI1zBQk9dVGExI/+KY+/0Dd49EkFohYUUd0zid3vC4wQAAk92MZ3e/xp5qIiIi8sBATV737o7jMFlcczhPG9MbI67yjQVcOkIhV+DR8b/HsDhXD7vJbsZzu9cyVBMREZEbAzV51an8Knx9MB8AEBqsxPzpQyWuqPNUChUeG/8H9lQTERFRixioyWtEUcRbnx5xv59100BoQ1QSVuQ9KoUKj47/vTtUm+0WPLv7NeSUMVQTEREFOgZq8pq9h4tw7FwlACAhJgS3jesrcUXe1RiqR+gvhurnvnsNh4uPS1wZERERSYmBmrzCanMg/Yuj7vcLf3Y1FD62vLg3qBQqPDLu9xihHwLAFapXf7cW3+Xuk7gyIiIikkrPSzwkic/2nEVJZT0AYMRV0UgdEidxRV1HpVDhkfG/x7XxwwEADtGJtfs24pNjX0IURYmrIyIiou7GQE2dVlVjxr/+exKAa0XEX8242r10d0+lkivxp3G/w039J7rbtvz0Kd7K2AKH0yFhZURERNTdGKip0zZ/mXNxmrzrktE3XidxRd1DJpPhV6Pvxezhd7rbdp7Zgxd/2ACz3SJdYURERNStGKipU3KLjNi5Lw8AEKxWYM4t/rciYmcIgoA7B9+MB66bD7lMDgDIKPwJT379d5TXV0pcHREREXUHBmrqlLc/OwJnw7DhX05NQURYkLQFSWRin+vw+MQHEKx0ff8zVXn4839W40jJCYkrIyIioq7GQE0ddvRsBbJOlgEAYiM1mDGhn8QVSWtY3CA8feOfEBMSBQAwWmrx9O5X8VnOTj6sSERE1IMxUFOHbd15sfd19k0DoVLKJazGN/QOT8Dz05ZjZMO0eqIo4v3sj/H3/70Fk80scXVERETUFfwyUOfk5GDmzJnQ6/WIjY1FamoqNm/e3ObPX7hwAc888wyuvfZaREVFITw8HAMGDMADDzyAoqKiFj8zefJkREZGQq/XN9vuuOMOb301v5GTV4nMht7puEgNJl2TKHFFviNUHYI/T1iMu4bc5m77Mf8QHv/v87hgbPnni4iIiPyX3wXqzMxMjBkzBg6HAzk5OSgpKcHixYsxb948rFq1qk3XGDJkCNasWYPnnnsOZWVlqKysxGuvvYYtW7Zg5MiROHOm5eWkP/74YxQXFzfbPv30Uy9+Q/+wdedJ9/7dU1J65CIunSGTyXDPsJ/h0fF/cI+rLjAW48//WY3/nN7NISBEREQ9iF+lIFEUsWDBAgDAxo0bER4eDkEQMH/+fMyZMwdPP/00srKyrngdp9OJ1atX46abboJMJoNMJsOtt96Kp556CqWlpVixYkUXfxP/djq/GgePlwAAYiKCceO1SRJX5LuuTRiOv05bjiRdPADA6rDhrYwP8Nc9r6PabJS4OiIiIvIGvwrUe/bsQXZ2Nm677TbodJ5zHc+aNQtOpxPr1q274nUee+wx/OxnP2vWPmHCBADA3r17vVNwD/VBk7HTM2+8CkqFX/0YdbteYbF4bupjuHnAJHdbZtER/OnLp3GwIFvCyoiIiMgb/CoJ7dq1CwAwevToZsca23bu3HnF6/zlL39BbGxss3ar1QoAiIqK6kyZPdq5QgP2HS0GAETpgjBtTG+JK/IPaoUKvxp9L/48YTF0QVoArllAXvh+Pf5xYDMfWCQiIvJjfhWojx8/DgBISEhodiwmJgZKpRJ5eXkwmUwdun5jz/SsWbNaPL5161aMHz8eSUlJ0Ov1mDRpEjZs2ACH4/JLTVssFhiNRo/NXzUdO33XDVdBqeDMHu1xTfzVePHmFbg2fri77b9nv8eyfz/F3moiIiI/5VeBurq6GgAQEhLS7JggCNBoNB7ntYfZbMa6deswaNAgPPDAAy2ek5ubiw0bNiA/Px/Z2dm44YYbsGjRIkyfPh02m63Va69evRo6nc69JSX555jjvGIjfjhcCAAID1PjpuuTJa7IP+mCtHhk/O/x22vnQC1XAQAqTFV44fv1ePGHDag0VUtbIBEREbWLXwXqrvToo4+irKwMH374oTuYN/Xhhx/iiy++wNChQwEAcXFxWLVqFWbPno0vv/wSb7zxRqvXXr58OQwGg3vLz8/vsu/Rlf7136a90wOg5rzTHSYIAqb2H4+XbvkLRjTMWQ0A+y9k4eF/P4n/nN4Np+iUsEIiIiJqK78K1OHh4QCAurq6ZsdEUUR9fb3HeW31/PPPIz09HTt27MDVV1/d4jkxMTGQy5sHyLvvvhsAsG3btlavr1arodVqPTZ/k19Sgz1ZBQAAXagKt1zfR9qCeojY0Gg8PvEBLLl+AbTqUACAyWbGWxkf4C+7XsSpinMSV0hERERX4leBevDgwQCAgoKCZsfKyspgs9mQnJyM4ODgNl/zb3/7G1avXo2vvvoK119/fbtrio93TYfW2oIwPcXH35xG49TJP580AEFqhbQF9SCCIGB88hj8/dYncEPfNHf7qYpzWPHfF/Daj+moqK+SsEIiIiK6HL8K1FOmTAEAZGRkNDvW2DZt2rQ2X++ZZ57B888/j127diEt7WKQOXjwoHvGDwDIysrCm2++2eI1CgtdY4pbmjWkpzDWWbE78wIAICRYiVvT+khbUA8Vpg7FH8bMxRM3PIyEML27fU/efjy04wn868gXMNstElZIRERELfGrQD1hwgSMGDECO3bsgMFg8Di2ZcsWyGQyLF682N0mimKr45VXrlyJ1157Dd98802zafhSU1PdQRlwBerHHnsMZnPzqc0++eQTAMCMGTM6/L183c59ebDZXeN5p6b2hiZIKXFFPdvQ2BT87ZaVWHjNPQhVuR7AtTps+Ojo/8PSHavwzdm9cDgvP7MMERERdR+/CtSCICA9Pd29YqLBYIAoiti4cSM2b96MlStXYuTIke7zlyxZgt69e+Ohhx7yuM4jjzyCZ599FjfeeCO2bduGVatWeWwtqaqqwty5c93DTWpra/H888/j3Xffxbhx47BkyZKu+tqScjhF7Phfrvv9beP6SFZLIFHI5LjlqslYc9uTuO2qGyAXXP9TrTRV440D72HZl0/h+7wDfHCRiIjIB/jdQNhRo0bhwIEDWLFiBVJSUuB0OpGcnIxNmzbhvvvu8zg3KSkJGo3GY5q66upqvPjiiwCADz74oE33vOuuuyAIArZt24aJEyeipqYGZrMZAwcOxAsvvIAHH3wQKpXKe1/Sh2QcL0Fppethz9GDYhEfHSpxRYElVB2C+df8EjcNmIh3sz/GocKfAABFNaVY8+M7+OT4l7jn6p8hNWEEBEGQuFoiIqLAJIhi46Nm1F2MRiN0Oh0MBoPPz/jxfxv2IvNkmWv/V9chdYj+Cp+grpRTdhof/PQZjpWd8mjvG5GEXwy5FakJIyAT/OofnoiIiPweA7UE/CVQF5TV4vd/dS33ro/SYP2fp0IuYy+o1ERRxJHSE/jgp8+aTasXHxaHGYOmYULyGCjlHOtORETUHRioJeAvgfrN7T/hsz1nAQALpg/FL24YIHFF1JQoisgsOoKtRz7HuSrPh28jgnW4PWUKpvYfD42y7dNIEhERUfsxUEvAHwK1yWLH/Ke+Qr3ZDpVCho1P3IwwTc8cJ+7vRFFEdvFxfJrzFY6WnvQ4FqwMwg19xuLmqyajV1jPndqRiIhISn73UCJ1j28PXUC92Q4AmHRNIsO0DxMEASN7DcHIXkNwquIcPj3+H+wvyALgWnVxx6lvsOPUNxjVayhuuWoyRuiHcJw1ERGRF7GHWgK+3kMtiiIefPEb5BXXAABeeXgS+ieGS1sUtUuBsRif5+zEnvMHYHPYPI71Co3F1P4TMKnPddAGhUlUIRERUc/BQC0BXw/UP50px+Ov/wAAGNwnEi88OEHiiqijaiy1+PrsXvzn9G6U1Vd6HJPL5Lg2fjim9BuH4XGDIZOx15qIiKgjGKgl4OuB+q/vHsAP2a6VIv80ZzQmXZMocUXUWU6nEwcLD+PLU9/iSOmJZsejNBGY3GcsJva5jmOtiYiI2omBWgK+HKgrDCYsfGYnnE4R4WFqvLPyJigV7LnsSQqNxfj63P+wO/dHGMzGZsf7RyZjQvIYpCWNRniwToIKiYiI/AsDtQR8OVB/uOsk3t1xHABwz9QU3HfrYIkroq5idzqQWXQEu87+gMyiI7j0rwJBEDA8bhDG9U7FtQnDEaoKkahSIiIi38ZZPsjD7kMX3PtTx/SWsBLqagqZHKkJI5CaMAKVpmr8kHcQe/L2Ibfa9TPQOB1fdvFxyAUZro4bhOsTRyE1cSS0ai5BT0RE1Ig91BLw1R7qc4UGLHnpWwB8GDGQXTAUYU/efnx//gDK6iqaHZcJMgyJuQpjEkfimvhhiA2JkqBKIiIi38FALQFfDdQbvziKbd+cBgD8/hfDcfu4vhJXRFISRREnys/ixwuHsO9CJirqq1o8L0kXj9Hxw3Bt/HAMiOzD2UKIiCjgMFBLwBcDtdMp4lfP7kR5tQkymYB3n7gZulC11GWRjxBFEWcq8/DjhUzsyz+EkrryFs/TqkMxXD8EI/VDMFw/GOFBvvHzTURE1JUYqCXgi4G66dzT1w6OwxO/vl7iishXiaKIvOoCHCr6CRkFh3G6Mg8iWv5rpE94Iobrh2CEfjAGRvWDSsEVN4mIqOdhoJaALwbqtR9m4asf8wAAf5wzGpM59zS1UbXZiMzCI8go/AmHS47DbLe0eJ5SpkBKdD8MjU3B1bEDMSCyDxRyPhdNRET+j4FaAr4WqG12B+au+gp1JhuCVHK8t+oWBKkZdKj97A47TlacbZgd5BjOVp1v9Vy1XIWU6L4YFD0Ag2IG4KqovghScJgRERH5H6YmwsHjpagz2QAA11/di2GaOkwhV2BIbAqGxKZg1vA7YDTX4HDJcRwpOYGjpSc9xl5bHFb8VHICP5W4Vm6UCzL0jeiNQdH9kRLdD1dF9UWUJkKqr0JERNRmTE6E3ZkX557mMuPkTdqgMIxPHoPxyWMAAGV1FThaetIdsCtMF2cOcYhOnK7MxenKXODkLgBAVHAErorq27D1QZ+IJPZiExGRz+GQDwn40pCPerMNc5/4Ela7E7pQFTb9382QyzntGXU9URRRXl+J42WnkVN+Bjllp3HBWHTZzwiCgERtL/SPSEb/SNfWOzwBKrmym6omIiJqjj3UAW7v4SJY7U4AwIQRCQzT1G0EQUBMSBRiQqIwsc91AIAaSy1OlJ/F6cpzOFVxDqcqcj0echRFEfmGQuQbCvFt7v8AuBaaSdDq0Tc8CX0iEtEnPBF9wpMQquZS6URE1D0YqANc06XGOdyDpBamDsW1CcNxbcJwAIDT6cQFYxFOVpzFmcrzOFOZi3xDIRyi0/0Zp+h0h+zv8va526OCI9A7PB5JugQk6xLQOzweCWF6zixCRERex/9nCWCVRjMOny4DAOijNBiYzAfAyLfIZDL0Dk9A7/AETO3varParcgzFOBMZR7OVOYht/oCLlwSsgGgwlSFClMVMouOutvkggz6sFgkanshSdcLiVrX1issFkoOGyEiog5ioA5g32UWwNkwgn7SqEQIgiBtQURtoFKo3A8qNrI5bMg3FCG3+gJyq/KRW52P84ZC1NtMHp91iE4UGItRYCzGvguZ7nZBEBAbEo2EsDjEh8UhXhuH+DA94sNioQvS8n8bRER0WQzUAWz3oXz3Pod7kD9TypXoF9kb/SJ7u9tEUUSFqQrnqwtw3lCI89UFyDcUoqCmBHan3ePzoiiipLYMJbVlOFR0xONYkEKNXqGx0IfFoldYDPShsdCHxiAuNAbhDNtERAQG6oBVUlmP0xcMAID+iTokxYVJXBGRdwmCgGhNJKI1kbgmfpi73eF0oLSuAvmGQlwwFuGCoQgFNcUorCmFpYVVHs12C85V5+NcdX6zYyq5EnEh0YgLjUFsaDRiQ6IaNtd+kDKoS78jERH5BgbqAHXwWLF7f+ywXhJWQtS95DI5eoXFoldYLMZgpLtdFEVUmqpRWFOCAmMxCo0lKKotRXFNKUrrK9DSDKNWhw35xiLktzLdX5gqBDEhUYgOiUSMJgrRmgjXe00kojURCFOHsoebiKgHYKAOUPuPl7j3UwfrJayEyDcIgoAoTQSiNBEYFjfI45jdYUdpXTmKastQXFOKkrpylNSWo6S2DKV1Fc2GkDSqsdahxlrX6hLsSrkSUcHhrvsGu+4dGRyOSE04IoPDERUcDm1QGGQCp7MkIvJlDNQByGyx46fTriWgo3RB6Bsv7eIyRL5OIVcgXqtHvLb5L59O0YlKUzVKaytQWleO0jrXa1ldBUprK1Bprm6xdxtwPUxZXFuG4tqyVu8tF2QID9IhIrjJFqRDRHA4woO0iAjWITxIC506DDIZgzcRkRQYqANQ9qky2BoWc0kdouc/ORN1gkyQucdqD8FVzY7bnQ5UmqpRXleJ8vpKlNVVoLy+ChX1laior0K5qQomm7nV6ztEp3sKwMsRBAFadRjCg7QIDwqDTq2FLigMuoawrQtqeK8Og1Ydyvm4iYi8iH+jBqADHsM94iSshKjnU8jk7ocVW1NvM6GivgoV9dWoNF2y1VehymyE0VwDES33dAOuMeAGsxEGsxF5bahLowx2h+swdSi06lBogxreq1xtYeoQ16sqBBplMH/5JiJqBQN1gBFFEQeOuQK1SiHD8KuiJa6IiDTKYGh0wUjSxbd6jt3pgMFsRJXJgCqzAVWmalSbjagyGVFtNqDaZES12bV/6SI3Lam3mVBvM6GotrRNNcoFGUJVIQhVhyBUFYIwVUiT9xrXfsNriEqDUJUGISoNNIpgDkUhoh6PgTrAnC0woNLo+ufl4VfFIEjFHwEif6CQyd0PTV6OKIqos9aj2mKEwVwDg9kVtI2WWhjNNTBYamA018BoqUW1xXjZ4SZNOUQnDBbX59tDgIBgZZArZCsbQrYqGCFKDUKUwdCoGl6VwdCoGl6VwdAog9z7cpm8XfckIupuTFMBpulwj2s53IOoxxEEwdVrrA5BovbKU2LaHXbUWOtgtNSgxlLrCt6WWtRYal2zlFhqUWNxvdY2zFpibmG+7taIEN294WWo6NB3UstVCG4SsIOVQa73iov7wYqmr2oEud+rEdTQrlaoOGMKEXUJBuoAc/AYx08T0UUKucI9e0hb2Rw21FnrGwJ3Heps9ahtfLXWodbieq2zmVBnrXdtNtdrW4ajXMrisMLisKLabGz3Zy+lVqgRrFAjWBGEIIUaQUo1ghRqqBWu10s3tVyFIKUaarkaQQoV1I1tCjXUDe9VciWDOlGAY6AOIFU1ZpzMd80U0KeXFrGRGokrIiJ/pJQrER6sQ3g7QjjgGo5icVhRbzU1BGwT6t2vrrZ6m9ndo11vrYep8b3d9Wq2WS77cOaVWOwWWOwWVKPz4bwplVwJtbwhcCtUDfsqqOSufVVjm1x5cV+hgkquhEquatiUl7Qpm2yu9xz+QuSbGKgDSMbxUjROh8vhHkTU3QRBcPf8RiK8Q9dwik6Y7RaYbRbU200w2czu0G22W9yvJpsZJrsZZpvF9Wq3wGwzw2S3uI43tDmcDq98N6vDBqvDhhprnVeu1xq5IIOyIWQrm4ZtmQIqhQpKmcLzuMz1qpQrXG0y137jq0quhELWyqtc4bqeTOHeV8gUnO2FqAUM1AHkYNPp8oYwUBOR/5EJMvdY6o6G8qbsDrsrbDeEbIvd6n5vbujNNtstsDisDe+trjaHtaG3u+HVYXO3W+1WmB2WVhf06QyH6ISjoSapKC4J2Y37je3Khn2FxzH5xfcy1/uL58kveW2pzfNVLpO3uK8Q5JDL5Az91O38MlDn5ORg5cqV+P777+F0OpGcnIylS5dizpw57brOe++9hzVr1iAvLw8ymQwTJkzAs88+i5SUlC69rxRsdicOnXBNjxWmUWJgcqTEFRERSU8hVyBUrkCoOsSr1xVFEQ6nwzX+uyFgW+02WB1WWB1WWByufYvd2tC77Xq12F3jxZu2WR022BxW1+edje8veXXavVr/5diddtiddqD7btluckHWELQbAndD0FbI5Je8bwzkMsgb2tyvTc6Ty2RN9puc476PHDJB1vB66fuL7XL3fWTNriEXZJA1vm/YVwhyyGSyhjbXLwr8ZcE3+V2gzszMxKRJkzBlyhTk5ORAp9Nh06ZNmDdvHk6dOoVVq1a16TorV67E6tWr8fbbb+P++++HwWDAggULkJqaij179mD48OFdcl+pHDtbAZPF9bff6EFxkMv4P0gioq4iCIKrZ1auQIiq659XcYpO2J0OWB1W2Bx2j6B9cb8xhLvaGo/ZnHZYHTbYnfZmx6xOG+wOO2wNIdrWZL+x3dZk396Nwf5yHKITDocTVodN6lK8rmnwljUJ302DuEwQXIG94RxZQ4CXCQLkDSHf/dmG403fN57b9FoyQQZ5w3kyQYBMkLv35U33Gz7TeI/mW0vHLm1r4X3jvSEgSBmE8CCt1H8UHgSxK/5NqouIoohRo0bh7NmzyM/Ph0538YGYefPmYfPmzcjIyMDIkSMve52MjAykpqZi7ty52LRpk7vdYDAgKSkJKSkpOHDggPu3QG/dt5HRaIROp4PBYIBW2z0/EG99egSffncGAPDIfaMxcVRit9yXiIgCR2PPvDuAN4RtV6+2w91mc9jhEB0X29znNJ7n8PhM477D49jFNte17LA7nbA77e7zHE4H7KLD89yGNofTAWcHZp0h6Y1JGIk/jf+d1GV48Kse6j179iA7Oxv33HOPR6gFgFmzZuG9997DunXr8Oabb172OuvWrYMoipg1a5ZHu06nw6233op//etf2Lt3L8aNG+fV+0rpwLFiAIBMJuCagbESV0NERD1R0555f+AUnXA6ne6AfTGAXwzqjiZB3vXeebFddLqPORv2HR6fd3q0X/ys5zUaa3A2XlN0Nuw3Xrdpm9P9y4D73KbXaDjPKYquezQc68zsOL7GF6ep9I+f+Aa7du0CAIwePbrZsca2nTt3dvo6//rXv7Bz5053oPbWfaVSUFaLwnLXk+eD+0QiVKOSuCIiIiLpyQQZZHIZFP4VhzrEKTrhFEV34HaFb9Edxt3hvUlAdzb+QiA22Xe3i3Be+plmx52tbC0dEy+576XnXHw/IKqP1P85m/Grn6Djx48DABISEpodi4mJgVKpRF5eHkwmE4KDg1u8Rn19Pc6fPw+VSoWYmJhmxxuvnZOT49X7SulAk8VcxnB2DyIiooDjGosMAHJArpS6nB7H9/rML6O6uhoAEBLS/GlsQRCg0Wg8zrvcNRrPvVTjtauqqrx2X4vFAqPR6LF1p8bhHgCQOkTfrfcmIiIi6un8KlD7q9WrV0On07m3pKSkbru3wymitt71lLM+SoPE2NBuuzcRERFRIPCrIR/h4eEAgLq65itRiaKI+vp6j/Mud43Gcy/VeO2IiAiv3Xf58uVYtmyZ+73RaOy2UC2XCXj1j5NRXFGH8moT568kIiIi8jK/CtSDBw8GABQUFDQ7VlZWBpvNhuTk5MuOY9ZoNOjduzfOnz+PsrKyZuOoG689aNAgr91XrVZDrVZf4dt1LX1UCPRR3l24gIiIiIj8bMjHlClTALjmkb5UY9u0adO8fh1v3ZeIiIiIeh6/CtQTJkzAiBEjsGPHDhgMBo9jW7ZsgUwmw+LFi91toigiPz+/2XUWL14MQRCwZcsWj3aDwYB///vfuOaaa5CWltbh+xIRERFR4PCrQC0IAtLT0yGKIhYsWACDwQBRFLFx40Zs3rwZK1eu9FitcMmSJejduzceeughj+uMHj0ajz/+ON5//31s2rQJoijCYDBg/vz5AIB33nnHY6xxe+9LRERERIHDrwI1AIwaNcq9LHhKSgpiY2Oxdu1abNq0CU8++aTHuUlJSdBoNC0+APjMM88gPT0da9asQVxcHFJSUiCXy3HgwAGMGDGiU/clIiIiosAhiKLYc9ai9BNGoxE6nQ4GgwFarVbqcoiIiIioE/yuh5qIiIiIyJcwUBMRERERdQIDNRERERFRJzBQExERERF1AgM1EREREVEnMFATEREREXUCAzURERERUScwUBMRERERdQIDNRERERFRJzBQExERERF1AgM1EREREVEnMFATEREREXWCQuoCApEoigAAo9EocSVEREREdDlhYWEQBOGy5zBQS6CmpgYAkJSUJHElRERERHQ5BoMBWq32sucIYmN3KXUbp9OJwsLCNv3G4y1GoxFJSUnIz8+/4g8F+R7++fk//hn6P/4Z+j/+Gfo3qf782EPto2QyGRITEyW5t1ar5V8ifox/fv6Pf4b+j3+G/o9/hv7NF//8+FAiEREREVEnMFATEREREXUCA3WAUKvVeOKJJ6BWq6UuhTqAf37+j3+G/o9/hv6Pf4b+zZf//PhQIhERERFRJ7CHmoiIiIioExioiYiIiIg6gYGaiIiIiKgTGKh7uJycHMycORN6vR6xsbFITU3F5s2bpS6L2unIkSNIS0uDIAjIzc2VuhxqA4vFgq1bt2L69OnQ6/WIiopCTEwMbr/9dvz3v/+VujxqA7vdjn/9619YsGABhg4divj4eMTExGDkyJF4/vnnUVtbK3WJ1AEPPPAABEHA/PnzpS6F2mjy5MmIjIyEXq9vtt1xxx1SlweAgbpHy8zMxJgxY+BwOJCTk4OSkhIsXrwY8+bNw6pVq6Quj9rAbDZjxYoVmDRpEk6dOiV1OdQOjzzyCO69914MHz4cp06dQkVFBTIyMmC1WjFt2jSsXbtW6hLpCsrLy3HPPfcgOzsbW7duRWFhIYqKirBkyRI8/vjjmDp1Kux2u9RlUjv897//xeuvvy51GdQBH3/8MYqLi5ttn376qdSlAWCg7rFEUcSCBQsAABs3bkR4eLj7N/I5c+bg6aefRlZWlrRF0hUtXboUR48eRXZ2NoYOHSp1OdQOTqcTEyZMwHPPPYewsDAAQO/evbFlyxYEBwfjT3/6E6qrq6Utktrk7bffxtVXXw0AUCgUWLhwIX7xi19g37592L17t8TVUVtVV1djwYIFuO+++6QuhXogBuoeas+ePcjOzsZtt90GnU7ncWzWrFlwOp1Yt26dRNVRWy1fvhzbt2+XbKl66rhbbrkFK1asaNYeHR2NQYMGwWKx4NChQxJURm0VFRWF77//HiNHjmx2LDk5GQBgMBi6uSrqqAceeAAjR47EwoULpS6FeiCF1AVQ19i1axcAYPTo0c2ONbbt3LmzW2ui9mv8P23yP9OnT2/1mNVqBeAKbOS7lEolxo0b1+Kx/fv3IygoCNddd103V0UdsW3bNnz11Vf46aefkJOTI3U51AOxh7qHOn78OAAgISGh2bGYmBgolUrk5eXBZDJ1d2lEAa28vBynTp3CkCFDMHz4cKnLoXZwOp04d+4cFi9ejEOHDiE9Pb3Fv2PJt5SUlOD3v/893njjDej1eqnLoQ7aunUrxo8fj6SkJOj1ekyaNAkbNmyAw+GQujQADNQ9VuPYzJCQkGbHBEGARqPxOI+IuseaNWtgt9uxZs0aCIIgdTnURl988QXCw8PRr18/fPnll/jnP/+Je++9V+qyqA1+/etf46abbsLMmTOlLoU6ITc3Fxs2bEB+fj6ys7Nxww03YNGiRZg+fTpsNpvU5XHIBxFRd/nxxx/x17/+FU899RSmTJkidTnUDtOnT4fRaERZWRk++OADzJo1CzfeeCO2bNmC0NBQqcujVrz99ts4dOgQjhw5InUp1AkffvghIiMjIZfLAQBxcXFYtWoVzpw5g/fffx9vvPEGlixZImmN7KHuocLDwwEAdXV1zY6Jooj6+nqP84ioax07dgzTp0/HkiVLWnxYkfxDTEwMHnzwQTz99NP44osvsHz5cqlLolbk5uZi2bJlePvttxERESF1OdQJMTEx7jDd1N133w3ANUZeagzUPdTgwYMBAAUFBc2OlZWVwWazITk5GcHBwd1dGlHAOXLkCG688UYsXLgQL774otTlkBc0PnS6fft2aQuhVu3YsQNOpxPz58/3WAjkF7/4BQDXmNzGNl+Zy5jaJz4+HgBQVFQkcSUM1D1W4z8nZ2RkNDvW2DZt2rRurYkoEGVmZuKGG27A73//e7zwwgvu9tzcXBQWFkpYGV3Jt99+iw8//LDFY43PoVRUVHRnSdQOixYtQk1NTbOFQD7++GMAwD333ONu85XV9qi5rKwsvPnmmy0ea/w7NDY2tjtLahEDdQ81YcIEjBgxAjt27Gg2T+qWLVsgk8mwePFiiaojCgwHDhzAlClT8MgjjzRbnXTVqlX4xz/+IU1h1CbffvstnnjiCTidzmbHGqcd5bR5RF0rKysLjz32GMxmc7Njn3zyCQBgxowZ3V1WMwzUPZQgCEhPT3evmGgwGCCKIjZu3IjNmzdj5cqVLS5WQETesXfvXkydOhVxcXGor6/HqlWrPDauVOofjh8/jt/85jfunjCr1YqPPvoIf/zjH6HVajmEh6gbVFVVYe7cue5hrLW1tXj++efx7rvvYty4cZI/kAgAgiiKotRFUNfJycnBihUr8P3338PpdCI5ORlLly7l0qt+Yv/+/e7fvCsrK2Gz2RAdHQ25XI45c+bgpZdekrhCas2dd955xXGZTzzxRLOea/Id5eXl+PDDD/HJJ5/g1KlTqKurg8lkQkJCAqZOnYo//elP6NOnj9RlUhs9++yzeO2112C1WlFVVYWgoCDodDrExsbi8OHDUpdHraipqcHHH3+Mbdu24ejRo6ipqYHZbMbAgQNx77334sEHH4RKpZK6TAZqIiIiIqLO4JAPIiIiIqJOYKAmIiIiIuoEBmoiIiIiok5goCYiIiIi6gQGaiIiIiKiTmCgJiIiIiLqBAZqIiIiIqJOYKAmIiIiIuoEBmoiIiIiok5goCYiIiIi6gQGaiIi8rBixQoIgoAffvihXZ87f/48BEHA66+/3kWVERH5JgZqIiLykJGRAZlMhpEjR7brc9u2bYNMJsMvfvGLrimMiMhHCaIoilIXQUREviM2NhaRkZHIyclp1+fGjRsHmUyGPXv2dFFlRES+iT3UREQEAFi6dCkEQUBZWRlOnDgBQRDc2/Hjxy/72cLCQvzvf//DzJkzu6laIiLfoZC6ACIi8g1jxozBPffcg61bt+KWW27BddddBwAQBAEpKSmX/ewnn3wCAG0a7uF0OvH666/j7bffRklJCQRBgF6vx7Rp07Bw4cIr3ouIyNcwUBMREQBg9uzZKCgowNatW/HAAw/g9ttvb/NnP/roI6SmpiIpKemK5/7f//0f/v73v2PHjh2YNGkSAGDnzp2YOXMmgoKCsGrVqo5+BSIiSXDIBxERuR06dAgAMGrUqDZ/pqysDHv27GnzcI+PP/4YAwcOdIdpAJg2bRqWLVuG6Ojo9hVMROQDGKiJiMgtMzMTsbGxiI+Pb/Nntm/fDofDgbvuuqtN5/fp0weZmZl44oknUFxc7G5/4okn8MADD7S7ZiIiqXGWDyIiAgDU1dVBq9Vi6tSp+Oqrr9r8uZtvvhmlpaXIzMxs0/lnz57F7NmzsW/fPshkMlx33XW48847MW/ePOj1+o6WT0QkGfZQExERACArKwtOp7Ndwz2qqqrwzTfftGt2j379+uHHH3/EwYMH8dhjj6G8vByPPfYYBgwYgM8++6wjpRMRSYqBmoiIAACHDx8GgHYt6PLZZ5/BZrO1ebgHADgcDgDA6NGj8dxzz+HkyZPYvn07bDYblixZ0q6aiYh8AQM1EREBACoqKgAAERERbf7MRx99hCFDhmDQoEFt/kz//v3x448/erTdcccdGDp0KCorK9t8HSIiX8FATUREAC7O7LFkyRIsX74cq1atwu7du1s9v6amxj3dXXs9/vjjyMvLAwCIoohPP/0UR48exf3339+x4omIJMR5qImICABw++23Y/Xq1fjHP/6BF198EXa7HQMHDmz1/C+++AIWi6Vdwz0A4PXXX8c///lPTJs2DbW1tRBFEYmJiXjllVfw29/+trNfg4io23GWDyIi6pC77roLP/30E06ePCl1KUREkmIPNRERdcj111/PIRpERGAPNRERERFRp/ChRCIiIiKiTmCgJiIiIiLqBAZqIiIiIqJOYKAmIiIiIuoEBmoiIiIiok5goCYiIiIi6gQGaiIiIiKiTmCgJiIiIiLqBAZqIiIiIqJOYKAmIiIiIuqE/w+oXSKIrZ0+lwAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 800x600 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "plt.plot(time,y_analytical(time),label=r'$y$')\n",
        "plt.plot(time,dy_dt(y_analytical(time)),label=r'$\\dot{y}$')\n",
        "plt.legend()\n",
        "plt.xlabel(r'$t$ / s')\n",
        "plt.ylabel(r'$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 the generalized $\\alpha$-scheme we find the approximation of the time-discrete rate as\n",
        "\n",
        "$$\\displaystyle \\frac{y_{n+1}- y_n}{t_{n+1} - t_n} = [(1-\\alpha) f(y_n) + \\alpha f(y_{n+1})]$$\n",
        "\n",
        "This can be re-cast in a residual form, i.e based on an approximation error:\n",
        "\n",
        "$$r(y_{n+1}) = \\displaystyle \\frac{y_{n+1}- y_n}{t_{n+1} - t_n} - [(1-\\alpha) f(y_n) + \\alpha f(y_{n+1})]$$\n",
        "\n",
        "where we would like to find the solution $y_{n+1}$ such that $r(y_{n+1}) = 0$.\n",
        "\n",
        "Due to the non-linear dependence of $r$ on the sought solution $y_{n+1}$, the equation cannot be solved directly. We thus use the concepts of a Taylor series expansion truncated after the first member to arrive at a recursive expression:\n",
        "\n",
        "$$r(y_{n+1}) \\approx r(y_n) + \\left. \\frac{\\partial r}{\\partial y} \\right|_n \\underbrace{(y_{n+1} - y_n)}_{\\displaystyle \\Delta y_{n+1}} \\overset{!}{=} 0$$\n",
        "\n",
        "Thus, with the help of the linearization we proceed from the known solution to the new one in a linear step\n",
        "\n",
        "$$\n",
        "    y_{n+1} = y_n - \\left( \\left. \\frac{\\partial r}{\\partial y} \\right|_n \\right)^{-1} r(y_n)\n",
        "$$\n",
        "\n",
        "Because this is only a linear step, we make an error (truncation error of the Taylor series expansion). Thus, we take the new guess for $y_{n+1}$ to repeat the above step recursively until we arrive at a satisfactory solution (linearization always around the newly found state):\n",
        "\n",
        "$$\n",
        "    y_{n+1}^{i+1} = y_{n+1}^{i} - \\left( \\left. \\frac{\\partial r}{\\partial y} \\right|_{i} \\right)^{-1} r(y_{n+1}^i)\n",
        "$$\n",
        "\n",
        "To perform this update, we need two incredients: a residual evaluation and a Jacobian evaluation (aka \"linear slope\")."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {},
      "outputs": [],
      "source": [
        "def residual(y_i,y_n,alpha,dt):\n",
        "    return (y_i - y_n)/dt - ((1.-alpha) * dy_dt(y_n) + alpha * dy_dt(y_i))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {},
      "outputs": [],
      "source": [
        "def df_dy(y):\n",
        "    return -2.*k*(y_eq - y)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {},
      "outputs": [],
      "source": [
        "def jacobian(y,alpha,dt):\n",
        "    return 1./dt - alpha * df_dy(y)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Note, that the Jacobian is constant for a given time step size (usually it is not!!). Thus, in the present case a distinction between Newton-Raphson, modified Newton, or initial stiffness methods is immaterial. In other words, an initial stiffness approach can be chosen here without any loss of convergence quality. We will demonstrate the full Newton-Raphson scheme below.\n",
        "\n",
        "The Newton update function is then as follows (returns $\\Delta y_{n+1}^{i+1} = y_{n+1}^{i+1} - y_{n+1}^i$):"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {},
      "outputs": [],
      "source": [
        "def Newton_Update(res,Jac):\n",
        "    return -res/Jac"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, for every time step we advance the system by applying the generalized midpoint Newton-Raphson update rule. Each new time step value then serves as the initial condition for the following linearization step.\n",
        "\n",
        "We will *recursively* apply the Newton update until the residual is below a specified tolerance:\n",
        "\n",
        "$\\left| r_{n+1}^{i+1}\\right| < \\epsilon_\\text{abs}$\n",
        "\n",
        "with a suitably chosen absolute tolerance set here to $10^{-6}$.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "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",
        "    #Jac = jacobian(y[-1],alpha,dt)\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",
        "        #Jac = jacobian(y[-1],alpha,dt)\n",
        "        while True:\n",
        "            #evaluate residual\n",
        "            res = residual(y_old,y[-1],alpha,dt)\n",
        "            #if residual is below tolerance, above maximum iterations, stop iterations\n",
        "            if (np.abs(res) < absolute_tolerance or i > max_iter): \n",
        "                break\n",
        "            #evaluate Jacobian\n",
        "            Jac = jacobian(y_old,alpha,dt)\n",
        "            #perform linear step\n",
        "            y_new = y_old + Newton_Update(res,Jac)\n",
        "            #update counter\n",
        "            i += 1\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": 24,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "application/vnd.jupyter.widget-view+json": {
              "model_id": "431b7397a3f14cfcb1106edef00f0067",
              "version_major": 2,
              "version_minor": 0
            },
            "text/plain": [
              "interactive(children=(BoundedFloatText(value=0.2, description=r'$\\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=r'$\\Delta t$ / s'),\n",
        "          alpha=widgets.RadioButtons(options=[0,0.5,1], value=1, description=r'$\\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=r'analytical')\n",
        "    ax[0][0].plot(time_n,y_n,label=r'numerical',marker='d',ls=':')\n",
        "    ax[0][0].legend()\n",
        "    ax[0][0].set_xlabel(r'$t$ / s')\n",
        "    ax[0][0].set_ylabel(r'$y$')\n",
        "    #error\n",
        "    ax[0][1].plot(time_n,y_n-y_analytical(time_n),ls=':',marker='d',label=r'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=r'relative')\n",
        "    ax[0][1].legend()\n",
        "    ax[0][1].set_xlabel(r'$t$ / s')\n",
        "    ax[0][1].set_ylabel(r'error')\n",
        "    #iterations\n",
        "    ax[1][0].plot(iter)\n",
        "    ax[1][0].set_xlabel(r'time step')\n",
        "    ax[1][0].set_ylabel(r'number of iterations')\n",
        "    \n",
        "    plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Tasks:\n",
        "\n",
        "- Compare the number of iterations with the Picard implementation (Exercise 3)\n",
        "- Investigate the influence of the time-step size\n",
        "- Repeat the analysis for $\\alpha = 1$ (Backward Euler)\n",
        "- How does the solution change for the modified NR method and the initial stiffness method? Discuss the result in terms of accuracy and iterations required to reach the result.\n",
        "- Repeat the analysis for a function with a non-constant Jacobian (one example is in lecture notes, but you may choose any other)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": []
    }
  ],
  "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.9.18"
    },
    "widgets": {
      "application/vnd.jupyter.widget-state+json": {
        "state": {},
        "version_major": 2,
        "version_minor": 0
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}