04b_Generalized_Midpoint_Newton_Nonlinear.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 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": 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": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgwAAAF9CAYAAABoEDD2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABK9UlEQVR4nO3deZhcZZ33//e3lq7qfUs66ewrCSRAgAYTBQmbKCLjOI46Oo74yOD4MI4wxuenIzjqoKMzbuNciuKo6CgKiiIoi+wIBkhC9n0hZO900vtW1VV1//6o6k6n0kl30tV1qro/r+uq69Q55z7nfKshXZ8+5z73MeccIiIiIqfi87oAERERyX0KDCIiIjIoBQYREREZlAKDiIiIDEqBQURERAalwCAiIiKDCnhdQC5761vf6h577DGvyxAREckWO9kKnWE4hSNHjnhdgoiISE5QYBAREZFBKTCIiIjIoBQYREREZFAKDCIiIjIoBQYREREZlG6rFBGRMSuRSHDkyBGam5uJx+NelzPiwuEwU6ZMIRgMnva2CgwiIjJm7du3DzNjxowZBINBzE46DEHec85x9OhR9u3bx8yZM097e12SEBGRMaujo4PJkydTUFAwqsMCgJlRXV1Nd3f3GW2fM4HBzK43s0fM7Ckze8nMHjWz84awXZmZ3WNmK8zsVTP7qpkF0trUmtnvzGx5qs2ykfskIiKST3y+nPkqHHHDCUW59FO6B/iZc+4q59xiYC3wlJlNGMJ2fufcxcBi4DLgi70rzcwHPAysc84tAa4APmZmN2f+I4iIiIxOuRQYnnfO3dtv/uvAOOAtJ9vAzBYCfwn8B4BzLgp8C7jVzEpSza4DFqX2h3OuBfg+cLuN9vNPIiIiGZIzgcE59660RV2paegUm10NdAMb+i1bARQCl/Zrs9M515zWZiow70zrPV1P7PgT96z+FUc7m7J1SBERkYzJmcAwgCUkw8BDp2gzC6h3zrl+yw71W9c7PZi2XXqbPmZ2s5mtNLOVDQ0Np1/1Sfzp9Zd5ZNvT1LfrgVYiIpJ/cjIwpC4V3AHc7pw7fIqmxUAkbVnvfNFptOnjnLvbOVfnnKsbP3786RV+CqWh5BWStmh7xvYpIiKSLTkZGIAvA687574+SLsOTrxk0TvfeRptRlxfYIgoMIiIyPBNmTKFb3zjG8ctW79+PeFwmE2bNmX8eDk3cJOZ3QqcDbx7CM13ATVmZv0uS0xMTXf2a/PWtO3S24y4slRgaFVgEBHJae/45O88Oe7DX/+L02q/ZMkSVqxYcdyyW2+9lZtuuolzzjknk6UBOXaGwcxuInlXw3udczEzm2VmV59ikydIdnBc0G9ZHckOky/2azPHzCrS2ux1zm3NWPGDKC1QYBARkcxJDwwPPvggq1ev5gtf+MKIHC9nzjCY2fuAzwI3Auem7ni8CKgFnky1uRN4J1DnnOt2zm00s98CnwI+ZGZB4BPAt5xzvd/MjwJrgNuAfzWzMuBm4M7sfLKkMl2SEBHJC6f7l75XFi9ezCc/+UkaGxspLi5m2bJlfO5zn6O6unpEjpczgQH4X5L1PJu2vH9UCpPsqNh//IQbgf82sxWAn2S4+FzvSudcwsxuAL5nZstT+7jbOXd3pj/AqagPg4iIZNJFF11EQUEBK1euZPXq1QQCAW655ZYRO17OBAbn3KCPznLOLQOWpS1rBT40yHYHgBuGVeAwHTvD0OFlGSIiMkqEQiEuuOACHn74YX7yk59w7733ntFTKIcqp/owjGaloWIAWiNtHlciIiKjxZIlS/jOd77DkiVLuP7660f0WAoMWdJ7SaI1qjMMIiKSGYsWLcLn851we+VIUGDIksJAmIAvQCQWIRqLel2OiIiMAj//+c/56Ec/yoIFCwZvPEw504dhtDMzSkPFNHW10BbtoDpQ4HVJIiKShxKJBA0NDdxzzz2sX7+e++67LyvH1RmGLCrTWAwiIjJMzz//PLW1tdxzzz088MADVFZWZuW4OsOQRbq1UkREhmvp0qUkEomsH1dnGLKoVMNDi4hInlJgyCKN9igiIvlKgSGL9IhrERHJVwoMWdT3xMpuBQYREckvCgxZ1Dfao84wiIhInlFgyKLeR1yrD4OIiOQbBYYsKguVAnoAlYiI5B8FhizSXRIiIpKvFBiyqKRfHwbnnMfViIiIDJ0CQxYV+IOEAyHiiThdPd1elyMiIqPAnj17KCkpYfXq1SN6HAWGLDv2mGtdlhARkeGbNGkSa9as4ZxzzhnR4ygwZJn6MYiISCYFAgHmzJlDKBQa0eMoMGSZAoOIiOQjBYYsK9UjrkVEJA8pMGSZnlgpIiKZMGXKFL7xjW8ct2z9+vWEw2E2bdqU8eMpMGRZ7/DQuiQhIiLDsWTJElasWHHcsltvvZWbbrppRDpABjK+RzmlY6M9KjCIiOSi99z3MU+Oe/977zqt9kuWLOG73/1u3/yDDz7I6tWruf/++zNdGqAzDFlXpksSIiKSAYsXL2bnzp00NjYSiURYtmwZn/vc56iurh6R4+kMQ5bpkoSISG473b/0vXLRRRdRUFDAypUrWb16NYFAgFtuuWXEjpdTgcHMCoAvAJ8C5jjndp+i7Y3Ap4FDaauWAO9zzv021a4ZWJPW5i7n3H0ZKfo0aeAmERHJhFAoxAUXXMDDDz/MT37yE+69916CweCIHS9nAoOZzQB+AWwD/EPc7CvOuXv67WMm8CrwWL82a5xzSzNT5fCV9T3iWk+sFBGR4VmyZAn/9V//xTXXXMP1118/osfKmcAAlAAfBKYAfzeE9g8OsOxG4JfOua7MlZVZJQXFGEZHtJN4Io7fN9RsJCIicrxFixbh8/lOuL1yJORMp0fn3Abn3I7TaN/snGvunTczIxk0fjQC5WWMz+ejpKAIh6M9qrMMIiJy5n7+85/z0Y9+lAULFoz4sXLpDMNwXQG0O+dWpC2faGb3AbVAlOSZibucc/Es19enNFRCW7SDtkgH5eEyr8oQEZE8lEgkaGho4J577mH9+vXcd192uuSNpsBwI/DDAZbvAG53zm1P9ZP4I3A2MGBXUjO7GbgZYNq0aSNSaGmoBNrqdWuliIictueff54rr7ySefPm8cADD1BZWZmV446KwGBmpcA7gH9OX+ecu77f+91m9lXgB2Z2h3OucYD2dwN3A9TV1bmRqLf3Tok23SkhIiKnaenSpSQSiawfN2f6MAzTe4CnnHNHhtB2J2DArJEt6eT0xEoREck3oyUw3MgAlyPM7CozuyFt8eTUdM9IF3UyGu1RRETyTd4EBjO708w2mFk4bfkcYAbJvgnppgKfMrOiVNty4DbgAefc4REu+aT0iGsREck3OdOHITXK4x+BitSiX5rZAefcu1LzYaCI5OWE/m4EfnqSux6eAi4CnjGzbqAUeJrkaJKe0fDQIiKSb3ImMDjnosDSU6xfBiwbYPntp9hmL/DxTNSXSXpipYhI7nDOkRzKZ/Rz7sz78ufNJYnR5FinRw3cJCLipWAwSFdXzg4OnHE9PT0EAmd2rkCBwQO9lyRaI20eVyIiMrbV1NSwf/9+Ojs7h/XXdz5IJBLU19dTXl5+RtvnzCWJseTYEyt1hkFExEtlZcnRdg8cOEBPT4/H1Yy84uJixo0bd0bbKjB4oDAQxu/zE4lFiMaiFAQKvC5JRGTMKisr6wsOcnK6JOEBMzvWj0FnGUREJA8oMHikTGMxiIhIHlFg8EiphocWEZE8osDgkd5LEi3dulNCRERynwKDRyoLKwBo6m72tA4REZGhUGDwSGVh8j7Yxs5mbwsREREZAgUGj1SlzjA0drV4W4iIiMgQKDB4pDcwNHU1e1qHiIjIUCgweKSq95KEAoOIiOQBBQaP9F2S6G4Z9eOXi4hI/lNg8EhBoIDigiLiibjGYhARkZynwOChqrAuS4iISH5QYPBQVVEFoDslREQk9ykweKiy79bKZk/rEBERGYwCg4d0p4SIiOQLBQYPHRuLQZckREQktykweEiXJEREJF8oMHioSoFBRETyhAKDhzQ8tIiI5AsFBg+Vh0rxmY/WSDs98R6vyxERETkpBQYP+Xw+KsJlADR1t3pcjYiIyMnlVGAwswIz+3czi5nZjCG032Jmz6a9PpHWptbMfmdmy83sVTNbNmIf4AzosoSIiOSDgNcF9EoFhF8A2wD/EDc75Jxbeop9+oCHgUedc3eYWTnwqpm1OufuHmbJGVGpsRhERCQP5NIZhhLgg8CPM7jP64BFwNcBnHMtwPeB283MMnicM9Z3p0Rns6d1iIiInErOBAbn3Abn3I4M7/ZqYKdzrrnfshXAVGBeho91RvouSXRr8CYREcldORMYzlCxmf3IzJ5P9V+43czC/dbPAg6mbXOo3zrP9V2S0BkGERHJYfkeGLYCdznn3gy8m+QliF/3W18MRNK26Z0vGmiHZnazma00s5UNDQ2ZrvcEGrxJRETyQV4HBufc3zrnVqTeHwHuAN5uZotSTTqAUNpmvfOdJ9nn3c65Oudc3fjx40eg6uP1PuJaz5MQEZFclteBYQA7U9M5qekuYGJam4lpbT1VFa4AkmcYnHPeFiMiInISeRsYzOxcM7spbfHk1HRPavoEMMfMKvq1qQP2Oue2jnCJQ1IYDBMKhIjEo3T2dHldjoiIyIDyJjCY2Z1mtqFfp8Zq4P+ZWXVqfQj4DPAKsCrV5lFgDXBbqk0ZcDNwZxZLPyUzoyqc7PioyxIiIpKrciYwpEZ5fBb4VmrRL83sN/2ahEl2VOwdP2EdyQ6Oj6S2e5HkHRDvcM7FAZxzCeAG4AIzWw48B9ydK4M29ertx6COjyIikqtyZqRH51wUWHqK9cuAZf3mG4F/GcJ+D5AMDTmrUndKiIhIjsuZMwxjWZWGhxYRkRynwJADNBaDiIjkOgWGHHDsiZXq9CgiIrlJgSEH6ImVIiKS6xQYcoAuSYiISK5TYMgBlalxGJq7W0kkEh5XIyIiciIFhhwQ8AcoC5XgnKM50up1OSIiIidQYMgRfZcl9JhrERHJQQoMOaLvTolu3SkhIiK5R4EhR/SO9ni0s8nbQkRERAagwJAjxhdXAdDQcdTjSkRERE6kwJAjJpSMA6C+44jHlYiIiJxIgSFHTCgeD0B9uwKDiIjkHgWGHNF7huFw+xGccx5XIyIicjwFhhxRUlBMYTBMV6ybtmiH1+WIiIgcR4EhR5gZE4pT/RjaGzyuRkRE5HgKDDlkQon6MYiISG5SYMghff0YdKeEiIjkGAWGHFKTuiRxSJckREQkxygw5JCJqUsSh3VJQkREcowCQw6p0eBNIiKSoxQYcsi4oip85qOxs5meeI/X5YiIiPRRYMghAZ+fcUWVOJyeKSEiIjlFgSHH6JkSIiKSixQYckyNnikhIiI5KOB1Af2ZWQHwBeBTwBzn3O5TtJ0AfAJ4MxADyoH7gf90zsX6tWsG1qRtfpdz7r5M1p4pfWcYFBhERCSH5ExgMLMZwC+AbYB/CJv8NXAtsNQ512ZmU4FXgULgc/3arXHOLc1stSNHlyRERCQX5dIliRLgg8CPh9j+MMmzCW0Azrm9JM8wfGBkysuO3udJaCwGERHJJTlzhsE5twHAzKYMsf39AyzuAkKZrCvb+p4n0ZF8zLWZeVyRiIhIbp1hyIQlJM8y9DfRzO4zs+fN7Ekz+0czG8olD08UFxRRXFBEJBahJdLmdTkiIiJADp1hGC4zuxqYBrw9bdUO4Hbn3PZUP4k/AmcDt5xkPzcDNwNMmzZtxOo9lQnF49gV3cPh9iNUhMs8qUFERKS/UXGGwcymA98FbnDONfdf55y73jm3PfV+N/BV4GNmVjXQvpxzdzvn6pxzdePHjx/Zwk+i97KEHkIlIiK5Iu8Dg5nVAL8D/t45t3oIm+wEDJg1ooUNgx5zLSIiuSavA4OZVQJ/AD7jnHsutezmfuuvMrMb0jabnJruyU6Vp6/3Mdcai0FERHJF3gQGM7vTzDaYWTg1XwI8CjwJNJhZnZnVAR/tt9lU4FNmVpTaphy4DXjAOXc4u59g6CbqDIOIiOSYnOn0mBrl8Y9ARWrRL83sgHPuXan5MFBE8nICJEd5fEPq9emT7PYp4CLgGTPrBkqBp0mOJpmzatSHQUREckzOBAbnXBRYeor1y4Bl/ea/BHxpkH3uBT6eoRKzprqwAr/5aOpqIRqLUhAo8LokEREZ4/LmksRY4vf5GV9cDcDhTj3mWkREvDeswGBmb8hUIXI8PYRKRERyyXDPMPwqI1XICY7dKaF+DCIi4r1B+zCY2UDPbIBk58MBBz+S4ZtYUgPAwbacvZlDRETGkKF0erya5FMk29OWG/DmjFckAEwpnwjAvtaDHlciIiIytMDwJ6C9d2Ck/sxsXeZLEoCpZZMA2NtywONKREREhhAYnHN/cYp112S2HOlVXVRJYSBMa6Sdlu5WyvUQKhER8dBpd3o0s8mDt5LhMjOmlNcCsLdFlyVERMRbZ3KXxKMZr0IGNLUsGRjUj0FERLx2JoHBBm8imTClXP0YREQkN5xJYHAZr0IGNE2BQUREcoSGhs5hfX0YWg/inHKaiIh4R4Ehh1WGyykOFtIR7aS5u9XrckREZAxTH4YcZmZM1WUJERHJAacdGJxz545EITIwdXwUEZFckJFLEmZWZmZ/aWYKExnWe2vlXt1aKSIiHjqjwGBmXzKzBjO72MyKgBXA/wLLzezvMlrhGNd7SWKfBm8SEREPnekZhiuAs51zK4APAJXADGAOcEtmShOAqX13ShzQnRIiIuKZMw0MXc65I6n37wN+7Jw74pw7BHRmpjQBKA+XURoqoaunm6NdTV6XIyIiY9SZBoZSM5tuZlcAlwP3AJiZHyjKUG2SMk2XJURExGNnGhi+BewAngR+5pzbbGaLgaeBDRmqTVKmpDo+7tGdEiIi4pFBH289EOfcvWb2DDDBObcmtXgP8DlgS4Zqk5Tefgw6wyAiIl45o8AA4Jw7CBzsN38A0J/AI6Bv8KZW/XhFRMQbGho6D0wtO9aHIeESHlcjIiJjkQJDHigJFVMRLiMSj3Kko9HrckREZAzKqcBgZgVm9u9mFjOzGUNoX2Zm95jZCjN71cy+amaBtDa1ZvY7M1uearNsxD7ACDp2WUL9GEREJPuGHBjM7Ftmdo2ZFYxEIamA8BwwCfAPcbN7AL9z7mJgMXAZ8MV++/QBDwPrnHNLSA449TEzuzlzlWdH3xDRulNCREQ8cDpnGDqA/wCOmtlDZvZRM5uawVpKgA8CPx5KYzNbCPxlqiacc1GSt3veamYlqWbXAYuAr6fatADfB243s7x66mbvGYbXm/d5XImIiIxFQw4MzrnPOucuAM4Cfge8BdhgZutTlwLenBq46Yw45zY453acxiZXA90cP+7DCqAQuLRfm53Ouea0NlOBeWdaqxdmVU0HYFfjHo8rERGRsehMHm990Dn3Q+fcXwHVwCdIXkL4HnDEzO43sxvNrDrDtaabBdS74x+wcKjfut5p+kX/9DZ5YWr5JIL+IAfbD9Me7fC6HBERGWOG1enRORdzzj3tnFvmnDsHuIBkP4S/Bj6ciQJPoRiIpC3rnS86jTbHMbObzWylma1saGjISKGZEPD5mVmRvAKkswwiIpJtGb1Lwjm32zn3Hefc251zX8vkvgfQAYTSlvXOd55Gm+M45+52ztU55+rGjx+fkUIzZU7qssSOxt3eFiIiImNOTt1WeZp2ATVpnRcnpqY7+7WZyPHS2+SNOdUzANhxdLendYiIyNiTz4HhCZIdHBf0W1YHdAEv9mszx8wq0trsdc5tzUaRmTS7agYA2xt3c3zXDRERkZGVN4HBzO40sw1mFgZwzm0Efgt8KrU+SLID5recc+2pzR4F1gC3pdqUATcDd2a3+syYWDKe4mAhLd2tHO1q8rocEREZQ84oMJjZ/EwXkhrl8VmSYykA/NLMftOvSZhkR8X+lyBuTG27AngZ+DPJJ2YC4JxLADcAF5jZcpIdMu92zt2d6fqzwcz6LkvsbHzd22JERGRMOdOnVa42s7uBzzvnMvKnbmrgpaWnWL8MWJa2rBX40CD7PUAyNIwKs6tmsPbQZnYc3c0bplzgdTkiIjJGnOkliUtI9h3YbmYfH86ATXJ6dKeEiIh4YdDAYGbnp/oH9HHOrXfOXQ3cBPwTsN7M3jZCNUo/vR0fdzXu0aOuRUQka4ZyhuEFkkMpn8A59yDJMw0/Idnn4JGR6N8gx1QWllNdVElXrJsDbfVelyMiImPEUALDOcCpetgVAatIhoZrgXVm9m0zK89AfTKAOamzDBqPQUREsmXQwOCc2+uci/dfZma3mtnPzWwbcJTkI6QvBv6L5GWKecAmM3vDCNQ85s1WPwYREcmyM71L4pPAcuAu4CVgVeouh14/NbP/D/gRxw+sJBnQe4Zh51HdWikikoucc8QTjmhPnGhPgmhPnEhPnJ5Y8n00dmx5NJagpyfe9z4ai9PTkzi2PJagJ9bbLrU+tZ+eWIKvf+LNFARH/t6DMwoMzrkB+zSk+THw5TPZv5zarKppGMbuln30xHsI+oODbyQiMsYlUl/gkZ44kWjaNPW+d3009epdH019QffNpwWBaL82va9ElgbkjcYSuRsYhqgBuHIE9z9mFQULmVw2kX2tB3m9eX/fYE4iIvksFk/QHY0TicbojsbpjsRS83G6ozEiPfHj1vd+eXdHYse+/I8LArHj5qOx7N5Z5vMZoaCPYMBPQdBPQcBHQdBPKOgnGPRREPBT0Df1E0ytL0htk5w/1i4YSO4jmNpXb/vCguyMbDBigcElH3bw3Ejtf6ybXTWdfa0H2dG4W4FBRLLKOZf6ok5+kXdFjr26I3G6Ij10pdb1fvF3pb78uyLJL/GuaKwvEPROY/GR/0IPFSS/sHu/uHvn+5YVHFve++Xdf31B33tfv/fHzwcDyW38/rx5+sKQjOQZBhlBc6pm8Nzul9h+9DXeOnep1+WISI5zzhGJxumMxOjs7qGzO0ZXd4zOSPJ9Z/exL/3O7p7jQkDvuu5+y0bidLvPIFQQIFzgJ1wQIFTgT3sfIBw69oXev01y2bH5gmDv8kDqyz/5JX78A47ldCgw5Kmzxs0CYEvDDo8rEZGR5pyjKxKjoytGR3cPHV09dHT30NnVQ0d3jI6unr4Q0NHdGwCOn+/q7snol3ww4KMwFCAcClAUCiTfF/gJp6aFvctCgeQXfWpZ/zbhUIDC3i/5UICCgE9f6DlMgSFPTS+fTHFBEQ2djRxuP0JNyTivSxKRU4jFE7R39tDeFaW9qyf5vjNKR1dPcr6rp+99/2lvGMjEl31B0E9ROPkFXxQOUBQOJr/Y+5YF+77oi8KBvveF/d4XpUJAYJSdbpfBZTQwmFm5c64lk/uUgfl8PhaMP4tX9q9hw+GtXKnAIJIVsXiCts4obR1R2jp7aO2I0t4Zpa0zmnyfCgNtndHktCu5visSH3znpxAu8FNcGKQoHKSkMEhROEBxYZDi8LH3ReEgxakgUDTAVF/yMhyZPsPQZGY7gM845x7I8L4lzcIJ85KBoX4rV856k9fliOSdRMLR0d1DS3uElvbkF35rRyQ17f9KLmvriNLRHTujY/kMigsLKClKfuGXFAYpKSpITZPzxYUFfeuKU6/eMKAve/FapgPDW4Czga+amd85d3+G9y/9LKyZB8CGw1txzunanwjQHYnR3B6huS1Cc3uElvbeaZSW1LLWjmjfNHGa5/p7v/jLioOUFYcoKQpSWlRAWXFB3/vkK0hJ4bFlhaEAPp/+jUr+OqPAYGbznXNb0pc7554EnjSzn3GSB1ZJ5kwum0h5uIzm7lb2tx1iSlmt1yWJjIhYPEFzW4TG1u6+aVNbhKa25Hxzv/fd0dM79V8cDlBWEqK8uICy4hDlJckv/7LiUGp67FVaXEBxOKgvfhmTzvQMw2ozuxv4vHOuKX1latkJyyWzzIwFNWfx5z0r2VC/VYFB8k484Whpj3C0pYujLd00tnbTmJoeTb1vauumtSOKG+KJgGDAR0VpiIqSUN+0vO9V0G8+GQqCAZ3qFxmKMw0MlwDfBLab2ReA76Y/oEqyY2HNPP68ZyUbD2/TeAySU+LxBI2tEY40dyVfLcnX0ebu1LSLxrbIkC4J+AwqSkNUloapLAtRVRY+br6iJERlWZiKkhBF4YAuz4mMgDN9lsR64Gozeyfwn8DHzOyTzrlHM1mcDG7hhGQ/ho2Ht5FwCXymv5YkOzq7e2ho6uJwUycNzV0cbkxOG5pSwaCle0hhoLykgKqyMNXlhalpmKqyMFXlYar6BYLRNmqeSL4ZVqdH59yDZvYIcBvwSzN7Efjngfo3yMiYUDyOcUVVHOlsZE/zfmZUquuIZEa0J059Y2fydbSDQ42dHG7q5HBjJ/WNXbR1RgfdR2VpiHEVhcde5WHGVRRSXZ6cryoLEQxkZxx8ERmeTNwlUQSsAn4C3AKsM7PvAXdoTIaRZ2YsrJnHs7uXs+HwVgUGOS3tXT0cOtLBwSMdHDjazqEjnRw82sGhox0cbek+5bbBgI+aykLGVxZRU1nE+MrC5HxF8n11eVhhQGQUOdO7JG4FLk69ZgNRYA3wX6npB4BNZvYu59zLmShUTm7hhFRgqN/K9fOu9rocyTHd0RgHGjrY39DOgYZ2DhxJvj94pIPWjpOfJfD5jJrKQiZUFTGhqjg1LWJCdRETKosoLwnpbgGRMeRMzzB8ElgO3AW8BKxyzvX/zfNTM/v/gB8BC4ZXogxmQc1ZAGxu2EE8Ecfv0191Y41zjqa2CPsOt7G3vp19h9vYd7idfYfbOdLcddLtQgV+aquLqR1XTG11MRPHFVNbXcTE6mLGVxSq34CI9DnTTo9DOe/9Y+DLZ7J/OT3VRZXUltRwsP0wu5r2MLd6ptclyQhxztHcFuH1Q63sOdTGnvq2vmlHV8+A2wT8Ru24YiaNK2Hy+BImjS9h0vhiJo0rpqosrDsKRGRIRvLhUw3AlSO4f+lnwYR5HGw/zIb6rQoMo0R3NMaeQ228dqCV3QdbeP1gG68faj3pZYSSwiBTJ5QypaaEKTWlTJlQwpSaEiZUFulMgYgM24gFBuecA5473e3M7AbgDqAL8AOfcM6tPEnbG4FPA4fSVi0B3uec+22qXTPJvhX93eWcu+9068tVC2vm8eTOP7Hh8Bb+8py3el2OnKbmtgi79rew60BLcrq/hYNH2gd8QmFxOMC0iWVMry1j6oQSpk8oY9rEUipKQzpbICIjJqceb21mFwH3Apc45zaZ2fXA42a2wDmXHgp6fcU5d0+/fcwEXgUe69dmjXNu6QiVnRMWTpiHmbGpYQedPV0UBQu9LklOorktwva9TezY18LOfc3s2Nc84B0Jfp8xfUIJM2rLmTGpjBm1ZUyfWMa4Cl1GEJHsy6nAAHwGeNw5twnAOfd7M6snebvmHQO0f3CAZTcCv3TOnbyn1yhUFiph/rjZbG7YwZqDG3njtDqvSxKSD0Lasa+ZbXua2Lanme17mzjcdOL/moUhPzMnlTNrcjmzJpUzc3I50yeW6rZEEckZuRYYrga+mrZsBXANAwQG51xz/3lL/tn1d8B7Rqi+nFY36Xw2N+xgxf61CgwecM5R39jJ5t2NbNndyNY9Tbx2oPWE0Q7DBX5mT6lg7tQKZk+pYM6UciaNK9EtiiKS03ImMJhZFVAOHExbdQh42xB3cwXQ7pxbkbZ8opndB9SSHDPiQZJ9GEbV8y8unnwe/7v2AV49uIFYPEbAnzP/eUeleMKx+0ALG3cdZdNrjWzefZTG1shxbXw+Y9bkcuZNq+SsaRXMnVbJlJpS/AoHIpJncukbpTg1jaQtj5AcTXIobgR+OMDyHcDtzrntZjYD+CNwNslLHccxs5uBmwGmTZs2xMPmhomlNUwtq2Vv60E2NWznvIlne13SqBJPOHbtb2b9jiOs33mUza8dpaM7dlyb0qICzp5RxfwZlcyfUcXcKRWEQ7n0z0xE5Mzk0m+yjtQ0lLY8BHQOtrGZlQLvAP45fZ1z7vp+73eb2VeBH5jZHc65xrS2dwN3A9TV1Q3xgbq54+Ip57N300Fe2b9GgWGYnHPsrW9jzbYG1u04woadR04ICDVVRSycVc2CWdWcPaOKKTUl6pAoIqNSzgQG51xj6vbHiWmrJgI7h7CL9wBPOeeODKHtTsCAWUDjIG3zSt2k8/nNpsdYuX8dH7nwffryOk2tHVFWbz3M6m2HWbOt4YS7FyZUFXHenHGcO2ccC2ZVU1M51JNfIiL5LWcCQ8qTQHpvvTrgN0PY9kYGGFnSzK4Cip1zD/VbPDk13XMGNea0WVXTqCqsoLGrmV1Ne5hdNd3rknJaIuHYub+ZlZsPs2pLPdv3NB039kFFSYhFZ43n/LnjOG/OeGqqFBBEZGzKtcDwFeBZMzvbObfZzK4j2VHxOwBmdifwTqDOOdf3p5+ZzQFmkOybkG4q8BEze9I512lm5SQfx/2Ac+7wiH4aD/jMR92k8/jjzudZsX+tAsMAIj1x1m5v4JWNh3hl4yGa2o51mwn4jXNnVXPhvBoumFfD9IlluntBRIQcCwzOuVVm9gGSD6/qHenx2n6DNoVJdoBM/w1+I/DTk9z18BRwEfCMmXUDpcDTwBdG4CPkhIunnM8fdz7Pyv3reN+5N3hdTk5o7+ph5aZD/Hn9QVZtOUy059j/KuMqCrn47AlcNL+G8+aOp1CdFEVETpBzvxlTlw4eOsm6ZcCyAZbffor97QU+nrEC88CC8WdRGAyzp2U/9e0NTCgZ73VJnmjvjLJ8/UFeWHeAddsbiMWPXWuYO7WCSxZM5JJzJjJzUpn6eoiIDCLnAoMMX8Af4ILahfx5z0pW7F/H9fOu8rqkrOns7uHljYf405r9rN56uC8k+AzOmzOOJefWsnhhLeMqNHS2iMjpUGAYpS6efF4qMKwd9YEhHk+wdvsRnl65l5c2HiQSTV5u8BmcP3ccly2azOKFtZSXpN+xKyIiQ6XAMEpdMHEhQV+ALQ07ONLZyLiiKq9Lyri99W08+coenlm197iOi2fPqOLyC6fwxvNqqSwNe1ihiMjoocAwShUVFHLR5PN4ae+r/Gn3K6PmkdfdkRgvrN3PH1/ew+bdx4bQqB1XzJV1U1l64RQmVhefYg8iInImFBhGsaUzFvPS3ld5fvfLvPPsa/O6Y9/e+jYeXb6bp1fs6RttsTDk59LzJ3PNJdOZP6Myrz+fiEiuU2AYxc6feA7l4TL2tx1iZ+PrzKme4XVJpyWecLyy8SC/f+E11u04NoDnvOmVXPuG6Vy6aLJugRQRyRL9th3F/D4/l027mN9ve4pndy/Pm8DQFYnx5Ct7eOhPOzl0NPkYkVCBn6UXTuFtS2Ywe0qFtwWKiIxBCgyj3OUzF/P7bU/x4p6VfGjRuwn6g16XdFJNrd387vmdPLZ8d99lhwlVRbzjsllcdfE0Sgpzt3YRkdFOgWGUm14xhRkVU9jdvI9VB9azeOqFXpd0gsONnTzwzHaeeGUPPbEEAOfMrOKdl8/mkgW1+DU0s4iI5xQYxoDLZyxm95pf89zul3IqMBw62sEvn9jKs6v2EU898WnJubX81RVzmDd99N0GKiKSzxQYxoBLp1/Mz9b+hjUHN9LS3Up5uMzTehqaurjvya08+coe4gmHz2csvWgK775yLtMnelubiIgMTIFhDCgPl7GodgGrDqznhddX8HaPRn5saY/wyye28tjy14nFE/gMrqybyvuumUftOI2dICKSyxQYxojLZyxm1YH1PLf7pawHhkhPnIee38mvntpOVySGGbz5gsn8zVvmMaWmNKu1iIjImVFgGCMumnQupQXF7G7ex7Yjuzhr3KwRP2Yi4Xhu9T5++shmjjR3JeuYX8ON1y9gRq0uPYiI5BMFhjEi6A9y1exLeXDz4zyy7ekRDwyvHWjhrgfW9Q3fPKO2jP/zjgVcMK9mRI8rIiIjQ4FhDLl2zuU8vOUJXtq3miMdjYwrzvydCB1dPfz88S384YVdJBxUlIT40NvP5oq6abo9UkQkjykwjCHVRZUsnnohL+5ZyWM7nuNvz//LjO7/xbUH+N5v19HcFsFn8I7LZvH+a+drwCURkVFAgWGMeftZV/HinpU8tfNPvHvBdYQDoWHvs6m1m7t+s47l6w8CycdLf+yvzmPmpPJh71tERHKDAsMYM6d6BvOqZ7H16C6ee+0lrp17+RnvyznHM6v28oMHN9De1UNhKMCH37GAa98wHZ8uP4iIjCoKDGPQdfOuZOufd/HI9qe5Zs5l+Mx32vtoaY/w3/ev4eWNh4Dk3Q+3vHsR4ysLM12uiIjkAAWGMeiSyYuoLqrkYNth1hzcyIWTzj2t7ddsO8w3f/Eqja0RisMB/v6d53Jl3VTMdFZBRGS0Ov0/LSXv+X1+3jZ3KQB/2Pb0kLfricX54UMbuOP7y2lsjXDOzCq+/ckruOriaQoLIiKjnALDGHXlrDcRCoRYX7+FXY2vD9q+oamLT3/nBR58bic+n/G3b53Pl//vpdRUFWWhWhER8ZoCwxhVUlDMNbMvA+C+Db8/Zdu12xu49ZvPsm1PM+MrC/nqP17Ke6+Zp3EVRETGEAWGMeyd899CKBBi9cENbDuy64T1zjkeeHo7n/v+n2ntiLLorPF889bLma9HT4uIjDk5FxjM7AYzW2Fmz5vZi2ZWN0j7LWb2bNrrE2ltas3sd2a23MxeNbNlI/sp8kNZuJS3n3UFAPdteOi4dZGeOP/5s1Xc84dNJBy85+qz+PzfL6G8ZPjjNoiISP7JqbskzOwi4F7gEufcJjO7HnjczBY45w6dZLNDzrmlp9inD3gYeNQ5d4eZlQOvmlmrc+7uTH+GfHP9vKt5bPtzrK/fyob6rSycMI+W9gh3/uhltrzeRGEowCfffyFvWFjrdakiIuKhXDvD8BngcefcJgDn3O+BeuCWYezzOmAR8PXUPluA7wO3m7r2U1JQzPXzrgbgvg0Ps7e+lWXffp4trzcxrqKQ//j4ZQoLIiKSc4HhamBl2rIVwDXD3OdO51xz2j6nAvOGsd9R47qzrqC0oJitR3byqR//hkNHO5kzpZyvf+LNegy1iIgAORQYzKwKKAcOpq06BJzqWczFZvajVJ+HZ83sdjML91s/6yT77F035hUFC7mw6o0AxMZv5pIFE/j3/3spVWXhQbYUEZGxImcCA1CcmkbSlkeAU93svxW4yzn3ZuDdJC9B/DptvwPtk4H2a2Y3m9lKM1vZ0NAw1Nrz2ovrDvDkI4aLFuAraWXplX7CoZzq3iIiIh7LpcDQkZqmd8MPAZ0n28g597fOuRWp90eAO4C3m9mifvsdaJ8MtF/n3N3OuTrnXN348eNP7xPkoadX7uE/frqCWMzHguLkWYafrXuArp5ujysTEZFckjOBwTnXCDQDE9NWTQR2nsauetvOSU13nWSf/duOSY8u3803f7GahIP3v2Ued7zz3cyumk5TVwsPbHrE6/JERCSH5ExgSHkSSB93oS61/ARmdq6Z3ZS2eHJquic1fQKYY2YVafvc65zbOrxy89dTK/bw3V+vBeDD1y/gb66dj9/n5yMXvg/D+MPWp9jXmt71Q0RExqpcCwxfAa41s7MBzOw6oBb4Tmr+TjPb0K9TYzXw/8ysOrU+RPLWzFeAVak2jwJrgNtSbcqAm4E7s/GBctGL6w7w7ftWA/CRGxbwrivm9K2bUz2DK2e9ibhL8ONX78M551WZIiKSQ3IqMDjnVgEfAH5qZs8DnwWu7TdoU5hkR8Xe8RPWkezg+IiZPQu8SPIOiHc45+KpfSaAG4ALzGw58Bxw91gdtGnVlnq+9rOVJBz8zVvm8c7L55zQ5m/O+wtKCopZX7+Vl/a96kGVIiKSa0x/QZ5cXV2dW7kyfViI/LVx11E+d/dyoj1x/uLNs/nIDQtO+ljqJ3b8iR+supfqwkq++bbPEQ7qFksRkTHgpAMa5tQZBhk5e+vb+LcfvkS0J85b3jD9lGEB4KpZb2JW5TSOdjXxs3W/zWKlIiKSixQYxoCW9ghf+J+X6OiOseTcWv7vu88/ZVgA8Pl8/MPFH8Tv8/PHHc+z5uCmLFUrIiK5SIFhlOuJxfnyPa9Q35gc7vmf338hft/QHqExo3IK7134DgDueuWntEc6BtlCRERGKwWGUcw5x3/fv4ZNrzUyrjzM7f/nDYQLTm8ExxvmXcO86lk0dbfwP6/+coQqFRGRXKfAMIr96qntPLNqH+ECP3d8ZDHV5YWnvQ+fz8cti28kFAjx5z0reXHPihGoVEREcp0Cwyj16pbD/O+jmzGDT37gImZNLj/jfU0sGc+HFv0VAP+z8hcc7WzKVJkiIpInFBhGoaMtXXz93uS4VR+4dj6LF9YOe59XzbqUC2sX0tHTxTf//D/E4rFh71NERPKHAsMoE48n+M+fraK1I8qis8bz11edlZH9mhn/95K/o7qokm1Hd3HP6l9lZL8iIpIfFBhGmZ8/voWNu45SVRbik++/CN8Q74gYirJwKcve9FGCvgB/3Pk8z+z6c8b2LSIiuU2BYRR5dcthfvXUdnwGy/62jorS9Kd6D9/squl85KK/AeB/Vv2CnY2vZ/wYIiKSexQYRomm1u6+fgvvf+t8zp09bsSOdeWsN3LN7MvoScT42ovfp7W7bcSOJSIiuUGBYRRwzvHdB9Ym+y3MHc9fX5mZfgun8uEL3sPc6pkc7Wziq3/6Lt2xyIgfU0REvKPAMAr8ac1+XtpwiMJQgI+/d1FG+y2cTMAfYNmbPsr4oiq2N+5O3jmRiI/4cUVExBsKDHmuuS3C936zHoD/844F1FQWZe3YlYXlfPbyj1NaUMzqgxv4/oqfoaefioiMTgoMee57v11HW2fyUsS1i6dn/fiTyiby6TffQshfwHO7X+IX63+X9RpERGTkKTDksRfXHuDFtQcIF/j5x/csGvQJlCNlbvVM/vlNf4/ffDy4+XEe2vJHT+oQEZGRo8CQp1o7otz1m7UA3Hj9AiZUZe9SxEAuqF3IP1z8QQB+tva3PLj5cU/rERGRzFJgyFM/f2wzLe1RFs6u5m1LZnhdDgCXz1zMP1z8QQzj3nUP8uuNf/C6JBERyRAFhjy0+2Arjy3fjc9n/MO7zsvKXRFDdeWsN3LLGz6EmXH/ht/zy/UPqSOkiMgooMCQZ5xz/ODB9SQcXPfGGUyfWOZ1SSd484w38E+LP4zPfPxm06P8ZM2vSbiE12WJiMgwKDDkmZc2HGTdjiOUFgV5/7XzvS7npN407WJuXfIR/ObjkW1P840//4BILOp1WSIicoYUGPJItCfODx/aCMAH3no2pUUFHld0aounXshnL/84RcFCXtm3hi8+801aulu9LktERM6AAkMeefC5ndQ3djJ9Yilv9WDMhTOxcMJ87rzqU30jQn72yf9gf+shr8sSEZHTpMCQJ462dPGrp7YB8Pd/cS5+f/78p5tSXsuXrv5/zK6azuGOo/zLE1/l5X2rvS5LREROQ/5864xx9z+5je5onMULJ3L+WeO9Lue0VRSW8/kr/pnFUy+kK9bN11+8m5+ueUDPnxARyRM5FxjM7AYzW2Fmz5vZi2ZWd4q2E8zsy2b2gpk9a2arzewzZhZIa9ecWt//9d6R/zSZ0dDUxR9f3oMZ/O3bzva6nDMWChRw25Kb+NCid+M3H7/f+iRffOabNHY1e12aiIgMIjB4k+wxs4uAe4FLnHObzOx64HEzW+CcG+jC918D1wJLnXNtZjYVeBUoBD7Xr90a59zSES5/xPz66W3E4gkuWzQ5J2+jPB1mxtvnXcXsqhl8c/kP2HJkJ596/Et8tO4DXDJlkdfliYjISeTaGYbPAI875zYBOOd+D9QDt5yk/WHgP51zban2e4H7gQ9kodas6H924X3XnOV1ORkzf/xs/uMt/8J5E86mLdLO1178Pt95+Sd0Rru8Lk1ERAaQa4HhamBl2rIVwDUDNXbO3e+c+2Xa4i4gNAK1eaLv7ML5k5mW52cX0pWHy/iXy/+RD1/wHoL+IM/tfollj9/JhvotXpcmIiJpciYwmFkVUA4cTFt1CJh1GrtaQvIsQ38Tzey+VL+IJ83sH83MP4xys6L/2YX3jqKzC/35zMfbzrqC/3jLvzC7ajpHOhv54rP/xXdf/imtkXavyxMRkZScCQxAcWoaSVseAYb0KEYzuxqYBnwxbdUO4Hbn3JuBm4B/Ar59kn3cbGYrzWxlQ0PDUGsfEb8axWcX0k0um8idV32K9y58B0FfgGd3L+fWRz7P07v+rGGlRURyQC4Fho7UNP1yQgjoHGxjM5sOfBe4wTnX3H+dc+5659z21PvdwFeBj6XOapDW9m7nXJ1zrm78eO9uX2xo6uKJl19P9l14yzzP6sgmv8/PXy24jq+99Q7OnTCP9mgH31vxv/zr099gZ+PrXpcnIjKm5UxgcM41As3AxLRVE4Gdp9rWzGqA3wF/75wbyohAOwHj9C51ZNXvnt9JLO647PzJTJ1Q6nU5WVVbWsPtl3+Cf1r8YcpDpWw9spPPPPEVvv3SjznS0eh1eSIiY1LOBIaUJ4H0cRfqUssHZGaVwB+Azzjnnkstu7nf+qvM7Ia0zSanpnuGXfEI6IrEeOKV5F/Uf3nFHI+r8YaZcen0S/jWdZ/nhvnXEPAFeOH1V/jEo5/n3nUP0h7tGHwnIiKSMbkWGL4CXGtmZwOY2XVALfCd1PydZrbBzMKp+RLgUZKBosHM6lIDPX203z6nAp8ys6LUNuXAbcADzrnDWfpcp+XplXvp7I5xzswq5kyp8LocTxUXFPG357+Lb73tX3nj1Ivoiffw4ObHueX3t3P/hofpiA56tUpERDIgpwZucs6tMrMPAD81sy7AD1zbb9CmMMkOkJaa/wTwhtTr0yfZ7VPARcAzZtYNlAJPA18YmU8xPImE4/cv7ALgHZfl7BWTrKspGcetb7yJtx+9il+uf4j19Vv49cZHeGTbM7z9rCt529wrKAkVD74jERE5I+ac87qGnFVXV+dWrkwfFmJkvbr1MP9693Kqy8P8z2evIZBHD5nKps0N27l/w+/ZeDj5QK5QIMSVM9/I2+ddRU1xtcfViYjkLTvZipw6wyD0nV247o0zFRZO4ezxc/nXK25j0+Ft/HbzY6w9tJlHtz/D4zueY8nUC3nb3CuYWz0Ts5P+vy8iIqdBgSGHHDjSzsrN9QQDPq5dPN3rcvLCOTVncU7NWexu2sfDW5/gz3tW8mLqNbNyKtfOWcqbptURChR4XaqISF7TJYlTyPYliR/8bj0PPb+Lqy6eyq3vuzBrxx1NjnQ28vj253h614u0pe6kKC4o4rLpl3DFzDcys3KqxxWKiOS0k56WVWA4hWwGhs7uHj78b3+kszvGN2+7fMzfHTFc0XgPy/es4vEdz7GjcXff8ukVU7hi5hLeNK2O8vDoHj1TROQMKDCciWwGhj+8+Brf+806zplZxVf/8bKsHHOs2N20l2deW86fXn+lb/wGn/k4d8J83jStjkumLKIoWOhxlSIiOUGB4UxkMzDc9q3n2LG3mWUfuIjLL5ySlWOONT3xHlYdWM+zry1n7aFNxFPPqAj6Apxfu4A3TF7ERZPPpaRAt2eKyJiluyRy2d76NnbsbaYoHGDxubVelzNqBf1BFk+9kMVTL6Qt0s7L+1bzwusr2Nywg5X717Jy/1r85mNBzTzqJp/HhZPO1S2aIiIpCgw54NlX9wHwpvMmEQrm/FO3R4XSUAlXz76Mq2dfRmNXMyv2reWV/avZeHg76+o3s65+Mz969T6mlNVy4aSFXFC7kHnVswj49U9GRMYm/fbzWCLheHbVXgCWXqRLEV6oKqzg2rmXc+3cy2mLtLPqwHpePbCBtYc2sa/1IPtaD/LQlicI+Qs4p+Yszpswn3MnzGdq+SSN8yAiY4YCg8c2vXaUw01djKsoZOGscV6XM+aVhkpYOnMJS2cuIRaPseXIjmR4qN/M3pYDrD64gdUHNwBQFirhnPFncU7NXM4ZP5cp5bX4TINticjopMDgsd7LEUsvnILPp79Wc0nAH2DhhPksnDAfgMauZtYf2sLa+s1srN9KU3cLL+17lZf2vQpAcbCQs8bNZv642cwbN5tZVdMIB0JefgQRkYxRYPBQtCfOC2v2A3CFLkfkvKrCCi6fuZjLZy7GOceh9gY2Ht7GxsNb2dKwk6NdTcedgfCZj6nlk5hTNYO51TOYVTmdKeW1BHzqpyIi+UeBwUMrNtXT0R1j9pRypk3UIEL5xMyoLa2htrSGq2dfCsCRjka2HNnJliM72H7kNV5v2c/rzft4vXkfT+16AUjewjmjYgozq6Yxo2IqMyqmMLV8koauFpGcp8DgoWd6OzteqOGKR4NxxVVcWlzFpdMvBiASi/Ja0x62H93Njsbd7GraQ317A9sbd7O93+iTZsakkglMrZjEtPLJTCufxLTySdQUj8PnU58IEckNCgweaWmPsHJzPT6Dyy+Y7HU5MgJCgQLmj5/D/PFz+pZ1RDt5rWkPu5r2sLs5eQZif+sh9rclXy/tfbWvbdAfZFLpBCaXTWRK2UQmlU5kUukEaktrdEZCRLJOgcEjL6w9QDzhuHB+DZVlYa/LkSwpLig6riMlJJ97sa/lAHtaDrC337Sxq7nvkka66qJKJpXWMLEk+aotHc+EkvHUFI9TmBCREaHA4JEX1qY6O2oY6DGvwB9kVtV0ZlUd/0jzzmhX8sxD6yH2tR7kQGs9B9rqqW9v4GhnE0c7m1hfv/WE/VWEy5hQPI7xJeOoKa5ifFE144uTr+qiSgr8wWx9NBEZRRQYPNDeGWXTa434fcbF50z0uhzJUUUFhcytnsnc6pnHLY8l4hzuOMKhtgYOtR/mUFsDB9sPc7jjCIc7jtLc3Upzdytbj+4acL/l4TLGFVVSXVTJuMJKqooqGVdUSVVhBVWFFVQWlhNUqBCRNAoMHli55TCJhOO8OeMoLtQvZjk9AZ+fSaUTmFQ64YR1iUSCxq5m6juO0NBxlMMdR/umRzobaexsoqW7lZbuVnY2vn7SY5SGSpLhIVxGRWE5leFyKgvLqQiXUREup6KwjIpwmcaZEBlDFBg8sGLjIQAuWaCzC5JZPp+PccVVjCuuGnB9IpGgqbuFI52NqQDRwtHORo52NdPY2URjVwtN3S20Rdppi7Rz8kiRFPIXUB4upTxUSnm4jLJwKWWhEspDpZSFSikLl1BaUEJZqITSUIn6V4jkMQWGLIvFE6zaUg/AJbocIVnm8/moTl2OmMfsAdskEglaI200djXT1N1KU1czTV0tNKUudbR0tdAcaaO5q4VIPMrh1BmMoQj6g5QVlFASKqa0oPjYtO9VRHFB0XHviwuKCPkL9NwOEY8pMGTZxl1H6eiOMXVCCbXjir0uR+QEPp+PisJyKgrLT9nOOUdXrJvW7jaau9toibTS2t1Oa6SNlkgbramzFP2nPfEejnY1cbSr6bRq8puPooIiSoJFFAULKSoopChYSHGw8Lj53ldhMExhIExRMExhal6hQ2R4FBiy7JVNqcsROrsgec7M+r6gJ5bWDNreOUckFqEt2kFbpIO2aDvt0Q7aI53JaTQ57eg/7emkM9pFJB7tu0wynHoLA8kgEQ6Gku+DIUKBMIWBEOFAiMJgmHDqfcifnIaDIUL+gmPLAwWEAiHC/gIK/AUaXEvGDAWGLHLOsWJj6nKE+i/IGGNmhINhwsEw44urT2vbnngPHT1ddEY76ezppqOnk45oF509nXT2dCVf0e7kNNZNV08XXT3J+a5YhO6ebiLxaF9bujL3uYL+YDI8BAoI+0MUBIKE/AWEAslAEfIXUOAPUhDo975vmnofCJ6wPOgPUuA79j7oD+o5JOIpBYYs2ne4nYNHOygrLmDe9IE7pYnIiYL+IBX+IBXhM3/mSjwRpyvWTXdPhK5YN1093XTHInTHIn3vu2LdRGLRY8tj3URikb5lkViE7ni0b1kkHqUn3kNPvAeiHRn8xAPzmS8VJAJ9ISLoCxD0ByjwBQn4AxT4gwR8AYK+AAF/gKAv2O99sm3v+t73va9gqk3A5yfgC6amve2S7/0+f9/7gC+A33y61DNG5FxgMLMbgDtI/g3gBz7hnFt5ivZlwLeBBan2TwCfdc7F+rWpBb4H1AAh4F7n3NdG7EOcxCupuyPqzp6AX4+yFskqv8/f17kyUxIuQTTeQzTekwwR8SiRWJRoPEok1kM0ngwayTbRVLto3/tovKcvdET7XlF64rHjptFEcppwieRxiGTsM2TCcUHCesOEP7Ws/3s/fku+96fa9i3vW+frW37CNr3rfX585k9NfX3TZJtjU1+q/bHtffh65wd477MT5+WYnAoMZnYRcC9wiXNuk5ldDzxuZgucc4dOstk9QIdz7mIzKwCeBb4I/Etqnz7gYeBR59wdZlYOvGpmrc65u0f4Ix1H/RdERhef+fr6NhAqGfHjxRNxovEeehIxYvEY0UQP0ViUWCJGTyKWCh6x1HwPPfFY8pXoIZaIpdYll8cS8eR+UvtKzifb9cRjxPvWJ5fHE3FiiXjfseKJePLlEsl9JGI5FmOGz7B+QcKXCiE+fD4/PjP81n96LGQct8wMf6q9L7W92bF99rbx9b0/2bLj561v3rh+3tVZGcE1pwID8BngcefcJgDn3O/NrB64heRZh+OY2ULgL4HzUu2jZvYt4B4z+7Jzrh24DlgEXJ1q02Jm3wduN7MfOOfcyH+s5MOmtuxuJOA3Lpg3PhuHFJFRxu/zU+jzU+h1If0kXIJYb3hIBYpYIk7Mpd7Hk9O4i/eFjvT3sXhqmogfv7/UsnjfNom+bXvDStwl+s0nt+9t17uvRCJBLDWf6LePhHPJ96l9JFL7S6T27XDJ/RKnx+sf9Cm8de7SMRkYrga+mrZsBXANAwSGVPtuYENa+0LgUuCxVJudzrnmtDZTgXnAlkwUPphVW+pJODhv9jiKwhrdUURGB5/5KPD7YBQOJ55IJI6FiH6hondZ7/L+QSThXL+27tj6E16uL9QkEuntBptPhpnedUFfdr7KcyYwmFkVUA4cTFt1CHjbSTabBdSnnSU41G9d73Sgffauy0pgeGWjBmsSEcknPp8PH77c+aL0WC716OjtiZR+GSwCFJ1im4Ha02+bobTpY2Y3m9lKM1vZ0NAwaNFD4ZyjrTMK6HZKERHJT7kUnHrvSUp/mk0I6DzFNgO1p982HUD6vVjpbfqkOkLeDVBXV5eR/g1mxpc+9iaa2rqpLA1nYpciIiJZlTNnGJxzjUAzkP4n+ERg50k22wXU2PE3Afduv7Nfm4H22b9NVigsiIhIvsqZwJDyJFCXtqwutXwgT5Ds4LggrX0X8GK/NnPMrCKtzV7n3NbhFiwiIjIW5Fpg+ApwrZmdDWBm1wG1wHdS83ea2QYzCwM45zYCvwU+lVofBD4BfCt1SyXAo8Aa4LZUmzLgZuDOLH0mERGRvJdLfRhwzq0ysw8APzWz3pEer+03aFOYZEfF/pcgbgT+28xWpNo/CXyu3z4TqdEjv2dmy1P7uDvbgzaJiIjkM8vSuEV5qa6uzq1cedJRqUVEREabkz63INcuSYiIiEgOUmAQERGRQSkwiIiIyKAUGERERGRQCgwiIiIyKAUGERERGZQCg4iIiAxKgUFEREQGpYGbTsHMGoDXM7jLccCRDO5vLNLPMDP0cxw+/QyHTz/D4cv0z/CIc+6tA61QYMgiM1vpnEt/uJacBv0MM0M/x+HTz3D49DMcvmz+DHVJQkRERAalwCAiIiKDUmDILj0hc/j0M8wM/RyHTz/D4dPPcPiy9jNUHwYREREZlM4wiIiIyKAUGLLEzG4wsxVm9ryZvWhm6hl8msyswMz+3cxiZjbD63ryiZldb2aPmNlTZvaSmT1qZud5XVc+MbM3m9lvzOyZ1L/j9Wb2Ca/rymdm9nEzc2a21Ota8omZ3WhmW8zs2bRX+UgeNzCSO5ckM7sIuBe4xDm3ycyuBx43swXOuUMel5cXUgHhF8A2wO9tNXnpHuCfnHP3ApjZV4CnzGyhc67e08ryx/uBNc65LwKY2SJglZntdM793tPK8pCZTQKWeV1HHvuKc+6ebB5QZxiy4zPA4865TQCpXy71wC2eVpVfSoAPAj/2upA89XxvWEj5OskBX97iUT356NvAN3tnnHNrgGZgjkf15Lv/Bv7d6yJk6BQYsuNqYGXashXANR7Ukpeccxucczu8riNfOefelbaoKzUNZbuWfOWc2+ScawMwM5+Z/T0QAX7lbWX5x8zeAfQAj3ldiwydAsMIM7MqoBw4mLbqEDAr+xWJALAE6AYe8rqQfGNmtwMHgNuA65xz+z0uKa+YWTHwJZI/Pzlz15vZ02b2gpn92swuHukDKjCMvOLUNJK2PAIUZbkWEczMgDuA251zh72uJ9845+4Eakl+6T1nZm/0uKR882/A95xz6X9EydDVAztIBtZLgd8Cy81s8UgeVIFh5HWkpumnfkNAZ5ZrEQH4MvC6c+7rXheSr1zSz4Hnga94XU++MLMLgDcA3/O6lnzmnHvUOfdp51x3av7nwHLg0yN5XN0lMcKcc41m1gxMTFs1EdiZ/YpkLDOzW4GzgXd7XEreMbMC51w0bfEm4CYv6slT1wOFwNPJE12EU8u/lfo9+VHn3FaPast3O4ERvV1fZxiy40lO/A9Zl1oukhVmdhNwHfBe51zMzGaZ2dVe15VHVqUu5/Q3CVAfhiFyzv2bc+5C59xS59xS4H2pVbemliksDEFqPJr0S9qTgT0jeVwFhuz4CnCtmZ0NYGbXkbwG+h1Pq5Ixw8zeB3yW5HX3c1MDh10DXOppYfmlFPh470xqfJV3Az/0rCIZq5YAH+mdMbPLgSuA747kQfUsiSwxsxtIdjTrIjnw0K3OuRXeVpU/zKwA+CNQAZwPvAwcGOB2QRmAmfUw8CXILzjnPp/lcvKSmb2f5OWHQiCemv4QuMvpF+lpM7NvAYtJ9mlYC2x3zv21p0XlCTN7G/CPJEOsj+S/7a855349osfV/+ciIiIyGF2SEBERkUEpMIiIiMigFBhERERkUAoMIiIiMigFBhERERmUAoOIiIgMSoFBREREBqXAICIiIoNSYBCRnGJmXzazb3tdh4gcT4FBRLLCzL5mZo8Noel1wCMjXY+InB4FBhHJlouBV07VwMwmAXOBZ7NRkIgMnZ4lISIjysyCQAcQ7Ld4k3NuwQBtbwLe6Zy7/hT7ez9wG9AOhIFnnHP/ktmqRSSdAoOIjCgz8wEXACtJPplwDxBxzjUN0PY3wBPOubtOsq9Jqe3Pcs7tMrPxwBbnXPWIfQARAXRJQkRGmHMuAdQCbcAK59yhk4SFIHA1p+6/MIHk4+FnpPbdQLLPg4iMMAUGEcmGC4C17tSnNN8M7HXOvX6KNmuA/wWeNLOnzexmYEPmyhSRk1FgEJFsWASsHqTNoHdHuKS/A84FVgFfAtaYWUUGahSRU1BgEJFsOB9YN0ib64A/nKqBmU02syXOuY3OuU8BC4BJwFWZKVNETkaBQUSyIQDMN7NJA50NMLOZJL/4XxxkP3OBr5pZIDXvAwzYkcFaRWQACgwikg2fBd4H7AP+fYD1byd5d0TPIPvZAuwElpvZs8DDwC3OubUZrFVEBqDbKkXEc2b2CPBr59yPvK5FRAamMwwikgueJXm2QERylM4wiIiIyKB0hkFEREQGpcAgIiIig1JgEBERkUEpMIiIiMigFBhERERkUAoMIiIiMigFBhERERnU/w/cqAJ5TDZTxQAAAABJRU5ErkJggg==\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 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": 7,
   "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": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def df_dy(y):\n",
    "    return -2.*k*(y_eq - y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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": 10,
   "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": 17,
   "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": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6c2ee68dc40a48cd8d8ccbb9d8e9cdca",
       "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",
    "- 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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}