{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "059e1855",
   "metadata": {},
   "source": [
    "# Formule de Simpson\n",
    "\n",
    "On rappelle la formule de Simpson pour le calcul approché d’une intégrale :\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "\\int_a^bf(x)\\mathrm dx \\approx \\Delta_x\\sum_{i=0}^{N-1}\\frac{f(x_i) + 4f(x_i+\\frac{\\Delta_x}2) + f(x_i + \\Delta_x)}6 \\equiv I(f;a,b,N) \\end{split}\n",
    "\\end{equation*}\n",
    "où $\\Delta_x =\\frac{b-a}N$ et $x_i = a + i\\Delta_x$.\n",
    "  \n",
    "1. Ecrivez une fonction `simpson_slow` qui calcule l’intégrale d'une function $f$ entre $a$ et $b$ avec $N$ pas avec la méthode de Simpson en utilisant une boucle (for loop).\n",
    "\n",
    "2. Ecrivez une autre fonction `simpson_fast` qui fait la même chose sans utiliser de boucle (on suposera que la fonction $f$ est vectorisée).\n",
    "\n",
    "3. Calculez l’intégrale de $f(x) = \\frac{1}{1+x^2}$ entre 0 et 1 pour $N=1000$ et comparez le temps entre les deux fonctions en ecrirant `%timeit` avant la commande.\n",
    "\n",
    "4. Calculez la valeur théorique $I^\\star$ de l'intégrale et tracez en échelle logarithmique le residu $|I^\\star - I(f,0,1,N)|$ par rapport à $N$. Qu'elle est la vitesse de convergence de cette intégrale ?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "501f2f54",
   "metadata": {
    "tags": [
     "sol"
    ]
   },
   "source": [
    "--------------------------------------------------------------------------------------\n",
    "\n",
    "## Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d72954aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1.\n",
    "\n",
    "def simpson_slow(f, a, b, N):\n",
    "    N = int(N)\n",
    "    dx = (b-a)/N\n",
    "    x = a + np.arange(N)*dx\n",
    "    integral = 0\n",
    "    for i in range(N):\n",
    "        integral += (f(x[i]) + 4*f(x[i]+dx/2) + f(x[i]+dx))*dx/6\n",
    "    return integral"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4bc63934",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2.\n",
    "\n",
    "def simpson_fast(f, a, b, N):\n",
    "    N = int(N)\n",
    "    dx = (b-a)/N\n",
    "    x = a + np.arange(N)*dx\n",
    "    return np.sum(f(x) + 4*f(x+dx/2) + f(x+dx))*dx/6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "29e3de20",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Temps d'exécution pour simpson_slow:\n",
      "5.59 ms ± 250 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "\n",
      "\n",
      "Temps d'exécution pour simpson_fast:\n",
      "60.6 µs ± 708 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "# 3.\n",
    "\n",
    "def f(x):\n",
    "    return 1/(1+x**2)\n",
    "\n",
    "print(\"Temps d'exécution pour simpson_slow:\")\n",
    "%timeit simpson_slow(f, 0, 1, 1000)\n",
    "print('\\n')\n",
    "print(\"Temps d'exécution pour simpson_fast:\")\n",
    "%timeit simpson_fast(f, 0, 1, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d874f58",
   "metadata": {},
   "source": [
    "La solution analytique de l'intégrale est : $\\int_0^1 \\frac{dx}{1+x^2} = \\arctan(x)|_0^1 = \\pi/4$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d09700e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "N_array = np.arange(6,30)\n",
    "residuum = []\n",
    "for N in N_array:\n",
    "    residuum.append(np.pi/4-simpson_fast(f, a=0, b=1, N=N))\n",
    "# Remarque : utiliser une boucle for ici ne ralentit pas le code car\n",
    "# 1/ le calcul de simpson prend relativement beaucoup de temps\n",
    "# 2/ il y a moins de 30 itérations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2d50c830",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fca3c2a9520>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEaCAYAAADpMdsXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABC0klEQVR4nO3dd3hUVfrA8e+bTkJooUgxQ5FeREEBpYko6NrbqlGxsu7Pvususti2sLrNtrpWkNUgrLoqGkRcC0UUBVEhVBEh9BIQEkJCyvv7496QYUjPTG4yeT/PM0/m3rn3nHfuzcyZc86954iqYowxxlRFhNcBGGOMqX+s8DDGGFNlVngYY4ypMis8jDHGVJkVHsYYY6rMCg9jjDFVZoWHaZBEZJ6I3FyN/TqKiIpIVCjiKifff4nInyrY5mERSQ1SfteLyGdBSitocQWkmyIiHwY7XVM5VnjUcSJytYgsFZFsEdkuInNEZKjXcZnaIyLjgTxVvd9v3UgR2eJhWLVCRIaKyOcisl9E9orIIhE5BUBVp6vq2V7H2FDV6q8nUzUi8ivgPuBWYC5wGBgLXAgE5VdhTYlIlKoWeB1HOFPVF7yOwQsi0gRIA34JvA7EAMOAPC/jMg6redRRItIU+ANwm6q+paoHVTVfVd9T1d+428SKyBMiss19PCEise5rI0Vki4j8WkR2ubWWG9zXBovIDhGJ9MvvYhFZ7j6PEJH7ROQHEckUkddFpIX7WnGzzU0ikgF8IiKRIvIPEdkjIj+KyO3+TTsi0lREprgxbBWRPxXnXdw8IiJ/F5F97v7n+MXVQkRedt/fPhF5x++180TkWxH5yf112q+c43mWiKxxf8E+DUjA6zeKyGo3j7ki4qvkebrB3S9LRDaIyC/K2fZ695fz427MG0TkNHf9Zvc8jfPbPtY9LhkislNEnhORRiKSAMwB2rk10mwRaefuFiMir7jxrBSRgX7p9RSnue4n97UL/F5LEpF3ReSAiHwFdAmIvYeI/M/99b9WRK4o5312EpH5bgz/A1oGvD7YPV8/ich3IjKyjKS6AajqDFUtVNVDqvqhqhb/nx7VtOb+z/2fiHzv5v1HEekiIl+47+t1EYlxty3+fPzO/b/dKCIpfmmdKyKr3HS2isi9fq/dIiLr3WPxrt+xL47hVjeGfSLyjIgc9b8WNlTVHnXwgVPDKACiytnmD8BioDXQCvgc+KP72kh3/z8A0cC5QA7Q3H39B+Asv7TeAO5zn9/tptsBiAWeB2a4r3UEFHgFSAAa4dSMVrnbNwc+creJcvd5x00jwY31K+AX7mvXA/nALUAkzq/MbYC4r88G/uOmGw2McNefDOwCBrn7jQM2ArGlHKeWwAHgMjeNe9xjc7P7+kXAeqAnTm38fuDzMo558fsvfm8/w/miFWCEe4xPLmPf6918b3Bj/hOQATzjHuezgSygsbv9E8C7QAsgEXgPeMTv/G4JSP9hINc915HAI8Bi97Vo9z3+DucX/Cg3r+7u6zNxft0nAH2ArcBn7msJwGY37ij32O8BepfxPr8AHnPf03A3n1T3tfZAphtjBHCWu9yqlHSauK/9GzgH93834Hh+5res7vFqAvTGqaF8DHQGmuL8j44L+HwUxzkCOOh3PLYDw9znzYvPqXvc9rjHIBb4J7AgIIY0oBmQDOwGxnr9fRKS7yivA7BHGScGUoAdFWzzA3Cu3/IYYKP7fCRwCL/CB+fLdrD7/E/AVPd5ovvB8bnLq4Ez/fZri/MFH0XJl2dnv9c/wS0M3OXR7jZRQBv3Q9zI7/WrgE/d59cD6/1ei3f3Pc7NtyjwS8Pd7lncgtJv3VrcwiVg/XW4X6LusgBbKCk85gA3+b0egVMI+EpJq/j9l1qo4xSUd5Xx2vXA937Lfd202vitywT6uzEeBLr4vTYE+NHv/JZWeHzkt9wLOOQ+HwbsACL8Xp/h7hPpnt8efq/9mZLC4+fAwoC8ngceKuU9JuN8KSf4rXuNksJjAvBqwD5zcb/US0mvJzDNPV8FOIVDG7/jGVh4nO63/DUwwW/5H8ATfscvMM7XgQfc5xnAL4AmAfFMAf7qt9zYPXYd/WIYGpDmfZX5zNe3hzVb1V2ZQEsp/6qedsAmv+VN7rojaejR/RE5OP/s4HygLxGnmesSYJmqFqflA952mxV+wilMCnEKgmKbA+LYXMZrPpxfvdv90nsepwZSbEfxE1XNcZ82Bo4H9qrqvmPfOj7g18VpuukeH/D+S41PnU91YIxP+qWzF+fLu30paR1FRM4RkcVuE8ZPOL+oW5azy06/54fceALXNcapScYDX/vF9YG7vjw7/J7nAHHu/1A7YLOqFvm9vgnnPbbCKeg3B7xWzAcMCjjWKTgFfKB2wD5VPVhOWpcHpDUU54fCMVR1taper6odcGpE7XBqZGUJPJalHdtipcVZ/P9zKc653OQ2wQ3xe39H3o+qZuN8Vv3/VwLPgX+eYcMKj7rrC5wmiIvK2WYbzoexWLK7rkKqugrnQ3AOcDVOYVJsM3COqjbze8Sp6lb/JPyeb8dpsip2fEBaeUBLv7SaqGrvSoS5GWghIs3KeG1yQIzxqjqjlG23+8fktkEHxviLgLQaqern5QXnFrz/Bf6O82u4GfA+Af0p1bQH58uut19MTVW1+IuoqsNhbwOOFxH/z3wyTvPUbpxf4ccHvFZsMzA/4Pg0VtVflpLPdqC52y9TVlqvBqSVoKqPVvQGVHUNTi2kT0XbVlJpcW5z81qiqhfi/Mh5B6cGAQGfOXf/JJzj2KBY4VFHqep+4EHgGRG5SETiRSTa/aX7V3ezGcD9ItJKRFq621flevrXgDtx2qXf8Fv/HDBZ3E5jN/0Ly0nndeAuEWnvftFP8Hsf24EPgX+ISBNxOuO7iMiIioJz950D/EtEmrvvf7j78ovArSIySBwJIvIzEUksJanZQG8RucT9FX4nR/9qfg6YKCK93ffbVEQuryg+nL6DWNwvX3E6+oNy6ahbQ3gReFxEWrtxtReRMe4mO4EkcS6sqIwvcZrBfusex5HA+cBMVS0E3gIedv/PeuH0IRVLA7qJyLXuvtEicoqI9Cwl7k3AUuD3IhIjzmXl5/ttkgqcLyJjxLnQIs7tvO4QmJY4nfS/Ln5NRI7HafJcXMn3XBnFcQ4DzgPecJdTRKSpqubj9JcVutu/BtwgIv3dHw9/Br5U1Y1BjKlesMKjDlPVx4Bf4XTg7sb51XY7zi8hcPotlgLLgRXAMnddZc3Aafv9RFX3+K1/Eqdt+UMRycL5sA4qJ50XcQqI5cA3OL++Cyj5wF2H80W7CtgHvEkZzRSluBanTXkNTp/N3QCquhSnk/1pN831OG3gx3Df2+XAozhNDF2BRX6vvw38BZgpIgeAdJwaWblUNQunIHrdjeFqnOMWLBNw3tdiN66PgO5u3mtwzt8Gt/mntOY6/1gPAxfgvK89wL+A69x0wPm/aozT5DINeNlv3yycQvFKnF/eO3COV2wZ2V2N8/+yF3gI5+KK4rQ241xq/jtK/qd/Q+nfRVluOl+KyEGc/8N04Nflvdcq2IFz3rYB04Fb/Y7HtcBG97jfClzjxv8x8ABOjXM7zsUSVwYpnnql+IoWY4LG/QX+nKr6KtzYGA+4Na9Uty/FVIPVPEyNiXPvwbkiEiUi7XF+bb7tdVzGmNCxwsMEgwC/x2kC+Abn6qwHPY3IGBNS1mxljDGmyqzmYYwxpsqs8DDGGFNlDWJU3ZYtW2rHjh2rte/BgwdJSEioeMMgq618Q5VPMNMNVlpenUsTXHYeS4T6WHz99dd7VLX0UQ28Hh+lNh4DBgzQ6vr000+rvW9N1Fa+oconmOkGKy2vzqUJLjuPJUJ9LIClamNbGWOMCRYrPIwxxlSZFR7GGGOqrEF0mBtjwkt+fj5btmwhNzfX61A81bRpU1avXl3jdOLi4ujQoQPR0dGV3qdeFh4ikowzIN4eYJ1WYjhnY0z42LJlC4mJiXTs2JFwneW1MrKyskhMLG0g6cpTVTIzM9myZQudOnWq9H613mwlIlPFmas5PWD9WHHmRl4vIvdVkEw3YLaq3ogzW5oxpgHJzc0lKSmpQRccwSIiJCUlVbkW50WfxzSc+bmPEJFInHmcz8EpDK4SkV4i0ldE0gIerXHGT7pSRD4BPg1VoKrKoq35HC4oqnhjY0ytsoIjeKpzLD0Z20pEOgJpqtrHXR4CPKyqY9zliQCq+kgZ+98LfKWqC0TkTVW9rJRtxgPjAdq0aTNg5syZVY5zzd5CHv0qlw6NhRv7xNK5WWSV06iu7OxsGjcO/eyVoconmOkGK63aOqYmtLKzs2nfvj0nnHCCp3E0a9aM3r17U1BQgM/n44UXXqBZs2ZVSmPZsmXMmDGDv/3tb8e81qdPH+bPn09SUlKZ+xcWFhIZGZzvpfXr17N///6j1p1xxhlfq+rAUnco6waQUD6AjkC63/JlwEt+y9cCT5ezfx+cCYWeA/5eUX41uUnwsf/8TwdN/kg73Zemf0pbqTl5BdVOqyrsJsHgp2U3l4WHTz/9VFetWlWlfVJTVX0+VRHnb2pqzeNISEg48vy6667TP/3pTzVP1I/P59Pdu3eXu82BAweCll9px5RybhKsKx3mpdWZyqwSqWo6ToFTfqIi5wPnt2/fnnnz5lUrsK7xuTx4SgKvr4vixYU/8s7SjdzYJ5aeSaGthWRnZ1c75rqQTzDTDVZatXVMTWhlZ2fTtGlTsrKyKrX9669HcccdcRw65HzNbNoEt9yi5ObmcsUVBTWKpTiGk046ifT0dLKystiwYQO//vWvyczMpFGjRvzzn/+kW7duvP322zz66KNERkbSpEkTPvjgAxYuXMhTTz3FG2+8QWZmJjfeeCOZmZkMGDCAoqIisrOz2bFjB1dccQVffvklAE899RTZ2dn87ne/45xzzmHy5MmcfPLJZGZmMmLECNLT05k+fTppaWkUFhayevVqbr/9dvLz85k5cyYxMTG8+eabtGjR4qj3kpubW7XPR1mlSigfHFvzGALM9VueCEwMVn7BGp7k8/V7dPhfP1HfhDS977/f6f5Dh6udblXyDSWreZj6pqo1D59PFY59+Hw1i6O45lFQUKCXXXaZzpkzR1VVR40apevWrVNV1cWLF+sZZ5yhqqp9+vTRLVu2qKrqvn37jryXn/3sZ6qqescdd+jvf/97VVVNS0tTQHfv3q0//vij9u7d+0i+f/vb3/Shhx5SVdWhQ4fqkiVLVFV19+7d6nPf1Msvv6xdunTRAwcO6K5du7RJkyb67LPPqqrq3XffrY8//vgx76e+1jyWAF1FpBOwFWdO4Ktrmmgwah6Bv1Z/dzK8sz6amV9tZs53WxjXO4aTWgf/MFrNI/hpWc0jPFS15pGR0ZjSGjcyMpSsrOxqx3Ho0CH69etHRkYG/fv3Z/DgwWzfvp3PP/+cSy+99Mh2eXl5ZGVlceqpp3Lttddy8cUXc/755xMZGUlOTg4FBQVkZWUxb948UlNTycrKYvjw4TRr1ozs7Gyys7MpKio68n7z8vKOpKmqHDx4kKysLLKzs1FVsrKyyM3NZejQoYBzD0eTJk0444wzyMrKomvXrkdqSf6qWvOo9cJDRGYAI4GWIrIFeEhVp4jI7cBcIBKYqqora5qXqr4HvDdw4MBbRo4cWa005s2bR+C+Y86E5Vt+4rdvLufJZVmcf2ISD5/fi6TGsTUNudx8QyFU+QQz3WClVVvH1ITWvHnziIuLq/T9DcnJTlPVseulRvdINGrUiOXLl7N//37OO+88XnnlFa6//nqaNWvG8uXLj9l+ypQpfPnll8yePZthw4bx7bffEh8fT1RUFImJiURERJCYmHgkJhGhcePGREU5X9PF61WV2NhYEhMTiY6OplGjRiQmJrJ//35EnPcUFxdH48aNj+wTGRlJUlISiYmJxMfHH8nLX1xcHCeddFKl33+tFx6qelUZ698H3g9mXqGoefi7t58yu3E07y7fxqertpHSM5bBbSODcgmh1TyCn5bVPMJDVWseDzxwdJ8HQKNGygMP5JKVVfM+j4iICB555BGuuuoqrrnmGpKTk3nllVe4+OKLUVXS09Pp27cvGzZsoFevXvTq1YtZs2axZs2ao2oegwcPZurUqfz2t7/lww8/ZN++fWRnZ9OkSRN27tzJxo0bady4MbNmzWL06NFkZWVx/PHHs2jRInr27ElqaupRNY/Dhw8fOUaqSnZ2NrGxsce8Vqxe9HnU9iPUQ7Kv3XFAL3z6M/VNSNMbXv5Kt/2UU+38qpJvMFifh6lv6uLVVqqq5513nr7yyiu6YcMGHTNmjPbr10979ux5pB/j4osv1j59+mjv3r31zjvv1KKioqP6PPbs2aNnnXWWnnTSSXr33XdrcnLykautnnzySe3SpYuOHj1ax40bd6TPY+nSpdq3b18dMmSITpo06ag+j9tuu+1IbP5XbgW+VqyqfR6ef7HXxqM25vMoKCzSlxZu0B73z9HeD36gqYs3amFhUcjzrSkrPEx9U53CI1x5eamuJzcJ1ha/ZqtbUlNTq5VGVW8s25VTxLSVeazKLKJ78whu6BPLcQlVv5HfbhIMflp2k2B4qCs3CdYFDe4mwdp+1PZMgkVFRTrzq03a56EPtNuk9/X5+es1v6Aw5PlWh9U8TH1jNY8SXtY8bD6PEBARfn5KMh/9agTDu7Xiz++v4ZJnP2f19gNeh2aMMUFhzVYVqGlTh6qyZEchr67OIycfzusczXldoomOKP+KLGu2Cn5a1mwVHqzZqoQ1W4VZs1Vp9mbn6d0zv1HfhDQd/Y95+vWmvbWSb0Ws2crUN9ZsVcKarRqA5gkxPP7z/rx8/Slk5xVw6bOf88e0VeQcrtl15sYY4wUrPGrZGT1a8+E9w0kZlMyUz35kzBMLWLR+j9dhGWNMlVifRwVC2U6+dm8hU9Pz2JmjDO8Qxc+7x5AQLSHP15/1eZj6xvo8SnjZ5xHWhUexgQMH6tKlS6u1b6jHQ8rNL+Txj9bx0sIfSUqI4awWfZj+l+PIyFCSk4XJkyElJWTZ29hWpt6ZN28ebdq0oWfPnl6HUmNFRUU88MADHDhwgIEDBzJu3Lgq7R+MOcyLrV69+phjKiJlFh7WbOWxuOhIJp7Tk3f+73Qi8mOZvulrDp64DGl0mE2bYPx4mD7d6yiNMaW58cYbad26NX369DnmtV/84hcsWrSI559/HhFh9erVR17r2bMnGzduZNasWWzdupXo6Gg6dOhQm6HXmBUedUTfDk3Z8e/T2begG/Fdd9Lu5vkk9N5CTo4yaZLX0RljSnP99dfzwQcflPral19+yeDBg1m+fDn9+/dn9uzZgDOk+s6dO/H5fKxdu5YhQ4bw2GOP8eyzz9Zm6DVmhUcdkrEpggNfdGXbtKHk702g5Xnf0fqyJWzdd8jr0IwxpRg+fPgxM/KB0wTUrVs3IiMjWbFiBffdd9+RwmPlypX07NkTEaFDhw40b94cIGh9F7XFCo86JDnZ+VuQmcjO6aex96NexB6/l/Y3z+fVLzZSVBT+/VPGhIM5c+YwduxYAFatWsUFF1zArl272L9/PytWrKBv374AXHLJJcydO5c77riD4cOHexlylYV1h3ldv9oq0Ecftebvf+9OXl7JL5D4Vln0u+krtmsu3ZpHcEPvWNo2Dl6Zb1dbmfom8Gqrv3z4A2t2Vn9GwNL0aNOYCWd3qdS2mzZtOmqOcYCLLrqIZ599lsLCQi644AKWLVvGhAkTGDRoEEuXLqVTp07ccsstNY7Ty6ut6so0tCGhIZpJMFRGjoSePWHSJPyutkrk6qtH8ebXW/hj2ioeWpzH3aO7Mn5YZ6Iia16I2NVWpr4JnEkwOiY66E0+0THRlb6KqXHjxkfNzJeTk0N2djbdunVj9uzZnHjiiSQmJnLRRRcxffp0tm/fzhVXXBGUq6SCebVVnZ9J0JQvJcV5zJs33++LTrh84PGM6NaKB2al89cP1jJ7+Xb+elk/erdr6mW4xnjuofN7ex3CUT799FPOOOMMgKOaqEaMGMGtt95KTk7OkXX1mfV51COtm8Tx/LUD+VfKyew8kMcFTy/ib3PXkJtf6HVoxjRIV111FUOGDGHt2rV06NCBKVOmHNXf4V94xMbG0rdvX2JiYmjWrJmHUQeH1TzqoXP7tuW0Lkn8MW01z3z6Ax+k7+Cvl/VjgO/Yqz6MMaEzY8aMY9adfPLJPP744wBMD7hJa9asWbUSV22wmkc91Sw+hn9ccSL/vvFUcvOLuOy5L3j43ZUczLOBFo3x0rJly4iOjvY6jJCzwqOeG9GtFXPvGc51g31M+3wjZz++gAXrdnsdljEmzNmluhXw6vLO6uS7bp8z0OKOg8rQ9lFc2T2GxjHeTDpll+qaULGBEUvYZFANYDKo2sr30OEC/cuc1dp54mwd8Mf/6fvLt4Ukn4rYZFAmVGwyqBI2GZQJmrjoSH47tgezbjud1omx/HL6Mn6Z+jW7snK9Ds2YoNIwbjWpbdU5llZ4hKk+7Zsy6/bT+c2Y7ny8Zhej/zGfN5Zutg+cCQtxcXFkZmba/3MQqCqZmZnExcVVaT+7VDeMRUdGcNsZJzC2z3Hc99/l/ObN5bz73Tb+fHFfFsyJd+9kH0FyMiGfN8SYYOrQoQNbtmxh9+6GfXFIbm5ulb/0SxMXF1flIeGt8GgAurRqzH/GDyH1y038Zc4azvz7AjI/7c7eTR0BOTJvCFgBYuqH6OhoOnXq5HUYnps3b16VhhQJJmu2aiAiIoTrhnRk7j3DObylBYkjVtEm5QuikrIAyMnB5g0xxlSaFR4NTIfm8WSknsKetBOJTsqm3fWf0WTweogoIiPD6+iMMfWFNVs1QMnJwqaVHTj0YytanLWS5iPWktBjOzHf9gNsoEVjTMXqZc1DRHqJyOsi8qyIXOZ1PPXN5MkQHw9FObHsmXUyu94aQGTjPCLGLOIvH9hAi8aYitV64SEiU0Vkl4ikB6wfKyJrRWS9iNxXQTLnAP9U1V8C14Us2DCVkgIvvAA+H4gorQ8fx0OnjOCyge15dt4PnPvkQpZs3Ot1mMaYOsyLmsc0YKz/ChGJBJ7BKRR6AVe5tYu+IpIW8GgNvApcKSJ/A5JqOf6wkJICGzfCJ5/MZ+NGuPm6aP562Ymk3jSIw4VFXP7cFzw4K51sG2jRGFMKT8a2EpGOQJqq9nGXhwAPq+oYd3kigKo+UkE6kcBbqnphKa+NB8YDtGnTZsDMmTOrFWt9GtsqWPnkFij//f4wH20qoEWcMK53DP1aVa17zMa2MqFi57FEqI9FnRvbCugIpPstXwa85Ld8LfB0Bfu/AEwHhlaUX0Ma2yqY+SzduFfP/Mc89U1I03v+843uzc4LSrpVZWNbGX92HkuE+lhQzthWdaXmcTkwRlVvdpevBU5V1TtqmE+DGlU3FPnkFynv/pDP+xvySYiGa3vFMrBNJCK1N1qv1TyMPzuPJazmAUOAuX7LE4GJwcrPah41z2fl1v163lML1TchTce/skR37j8UlHQrw2oexp+dxxJW8xCJAtYBZwJbgSXA1aq6sob5WM0jiPkUFilzN+bz9vp8oiPgyh4xDGsfVWotxGoeJlTsPJZoUDUPYAawHcgHtgA3uevPxSlAfgAmBTNPq3kEN58fdmXp5c99rr4JaZry4mLNyDwYlHTLYjUP48/OYwkvax61foe5ql5Vxvr3gfdrORxTDZ1bNWbmLYN57asMHp2zhrMfX8BvxnRn3GkdiYwovy/EGBMebBraClizVfkyDxXx71WHWb67kC5NI+iZ3ZY3p/Rg165YWrfO4+abNzB69C5PYwx2OsZbdh5LNKhmKy8e1mwV2nyKior07WVbtOfv5mryve9r09PWKRGFCqrx8aqpqd7HGMx0jLfsPJZocB3mtcVqHrWbz5XXn0xhvx9I6LWNw7sSyZzTj8M7mtGmTS4zZy72PEb7xRoe7DyW8LLmEdaFR7GBAwfq0qVLq7XvvHnzGDlyZHADqkP5BjOfiAhQhUYn7KTF2SuITMjjwJLOHFjUjcLDkZ7H6NW5NMFl57FEqI+FiJRZeNTLUXVN3ZSc7Pw9tL4N26aMIHv58TQdtIHjb1nA4g2Z3gZnjAkqKzxM0BQP9Q6gedHsnduPn94eRIskuPKFxUx6ewVZufneBmmMCYqwbrayPo/az+ejj1rz0kudj7raatjInby1/jAfbiygeZxwXa8Y+reu/FXi1udh/Nl5LGFXW9nVVp7lG6p8Skt32aa9etZjzkCLd81YppmVHGjRrrYy/uw8lvDyaitrtjK15qTk5qTdMYy7zuzK7BXbGf3YfN79bhsaxrVfY8KVFR6mVsVERXDPWd1Iu2MYxzdvxJ0zvuGWV75mx/5cr0MzxlSB9XlUwPo8QpdukSofbizgre8PExkBP+8ew4gOxw60aH0exp+dxxLW52F9Hp7lW5t9HmXZuCdbr3z+C/VNSNMrn/9CN+7JrnZawYrJ1F12HktYn4dp0HxJCbx2yyAeuaQv6Vv3M+aJBby0cAOvpiodO8KoUSPo2BGmT/c6UmNMsVofVdeY0ogIV52azBndW3P/Oyv40+zV5O/Yzu6D/VBNZNMmGD/e2TYlxdtYjTHWYW7qmOOaxvHidQPh85OIaJJD2+sX0vT0dRBRRE4OTJrkdYTGGLAO8wpZh7k36Y4aNQKJy6f5mStp3Hsbh3c7Ay3m72jKJ5/M9yQmUzfYeSxhAyOGmA2MWPv51DTdjh1h0ybneaPOO2kxJp3IxrnIuk6sfr07jWKqPtCiDagXHuw8lrCBEY0J4D9O1qENbdg2ZTi5K5Ohx4+MeWIBn/+wx9sAjWngrPAwdVJKCrzwAvh8IKIkt43msav7MnP8YCIErn7xSya+tZwDNtCiMZ6wwsPUWSkpsHEjfPLJfDZudJYHd05izl3D+cXwzvxnyWbOemw+H63a6XWoxjQ4VniYeqdRTCQTz+3JO7edTvP4GG5+ZSl3zPiGzOw8r0MzpsGwwsPUW/06NOPd24fyq7O68UG6M9DirG+32kCLxtSCsL7ayi7V9S6fYKZbmbS2ZhUxNT2PH/YXcWKrSK7rFUNSo6N/G9klnuHBzmMJG9vKxrbyLN+6MLZVsNIqKCzSKQs3aI/752jvBz/Q1MUbtbCwKCQxGe/YeSxhY1sZEwSREcKNQzsx9+7hnHh8Uya9nc5VLy7myakHbYwsY4LMCg8TdpKT4km9aRB/ubQv3206wGOrFrCvzQ8oemSMLCtAjKkZKzxMWBIRfn5KMrlvjeDQhlY0P2MNx137OdGtDtgYWcYEgRUeJqxlrItj99sD2P3OyUQ1OUTbcZ/RdOg6MrYWeh2aMfWaFR4mrCUnAwg5a9uy7aURHFzdjmanf0/yzZ+xLGOf1+EZU29Z4WHCmv8YWUW5MWTO7s/+d0+hacsCLn32c/7w3ipyDhd4G6Qx9ZAVHiasBY6R5fPB05Na89mk4VwzyMfURc5Ai599bwMtGlMVdb7wEJHOIjJFRN70W5cgIv8WkRdFxOaVM+UqbYysxLho/nhRH/4zfjBRERFcM+VLfvvmd+zPsYEWjamMkBYeIjJVRHaJSHrA+rEislZE1ovIfeWloaobVPWmgNWXAG+q6i3ABUEO2zQggzonMeeuYfxyZBf+u2wrox+fzwfpO7wOy5g6L9Q1j2nAWP8VIhIJPAOcA/QCrhKRXiLSV0TSAh6ty0i3A7DZfW6XzZgaiYuOZMLYHsy67XRaNo7l1tSvuW36MnZn2UCLxpSlymNbiUgE0FhVD1Ry+45Amqr2cZeHAA+r6hh3eSKAqj5SQTpvqupl7vNrgX2qmiYiM1X1ylK2Hw+MB2jTps2AmTNnVvYtHsXGtvI+3WClVZl0CoqUOT/mM2t9PrFRcGJhSz6Z2o9du+Jo3TqPm2/ewOjRu2oci6k+G9uqRJ0f2wp4DWgCJABrgO3Abyq5b0cg3W/5MuAlv+VrgafL2T8JeA74AZjorksAXgaeBVIqisHGtqr9fLwY2yqY6Xy/M0uH/X6R+iakaevLv9TIxBwF1fh41dTUoIRjqsnGtirh5dhWlap5iMi3qtrf7ZweAEwAvlbVfpXYtyNH1zwuB8ao6s3u8rXAqap6R4WBVJGNqutdPvW15uHvyisHkdN+B81GrAWFffN7kP2NjzZt8pg5c3GN4zHVYzWPEvWh5rESiAbeAEa4676r5L4dObrmMQSY67c8EbdGEaqH1TxqP5/6XvNQVRVRBdXIJge19RWL1TchTdtcvUijk7KCEo+pHqt5lKgPNY87cWob3wE/A5KBVFUdVol9O3J0zSMKWAecCWwFlgBXq+rKCgOpIqt5eJdPeNQ8BrNzZ5y7pCT02UrzUauIiCnk8p6RnNMxmsgIqXFcpmqs5lGiztc8SnsAUZXYZgZO/0g+sAW4yV1/Lk4B8gMwqboxVPZhNY/azyccah6pqU4fB5Q8Eloe0nMnL1XfhDT92VMLNH3rT0GJzVSe1TxKeFnziKpM6SMiD5bx0h/K209Vrypj/fvA+5XJ2xivpLi3n06aBBkZzjhZkyfHkZIygDkrtvPArJVc8PQibh3RmTtGdSUuOtLbgI2pRZVttvq132IccB6wWlVvDFVgwWDNVt7lEw7NVhWmd1iZufYwn20toG2CcGOfWLo2twIk1KzZqkS9a7YCYvHr9K7rD2u2qv18wqHZqrLmr92lpz3ysXa8L00fmpWu2bn5IcnHOKzZqkR9nIY2HuhczX2NCSvDu7Xiw3uGM25IR/79xUbOfnwBC9btZvp06NgRIiKw6W9N2Klss9UKoHjDSKAV8AdVfTqEsdWYNVt5l09DaLYqzff7CpmSnseOg0rOyvZkftSLotwYAGJjC7n33rV2h3oNWbNViTrfbAX4/B7tqcSVVnXpYc1WtZ9PQ2q2CnTocIH6zlutyb+ZrR1u+5/Gd9t25Gotn69WQghr1mxVos42W4lICxFpAWT5PQ4BTdz1xpgAcdGRZMzuwfZ/n05BdiytLl5Gy4u+JiIhl4wMr6MzJjjKbbYSkR9xmqsE58bAfe7zZkCGqnaqhRirzZqtvMunoTZbFTtyg2FEEU1O2UCzod9TVBBB0dKupP5xEyJ2c2F1WbNVifrQbPUccK7f8jnAPyqzb114WLNV7efTkJutVI+9wTCqRZa2u9YZaPGalxZrRubBWosl3FizVYk622zl5xR1buwrLnDmACNqUqIZE86Onv4W2ic25tGzh/DHC3uzbNM+xjyxgGmLfqSoqGpTIhhTV1TqDnNgj4jcD6TiNGNdA2SGLCpjwkBKSsld6g4BOnJGj9ZMejudh99bxXvLt/OXS/tyQutEj6I0pnoqe6luC+AhYLi7agHwe1XdG8LYasz6PLzLp6H3eVREVfl8WwGvrTlMXgFceEI053SKJsoGWqxQXTqPXqvzfR71/WF9HrWfT0Pv86isXQdy9f9Sv1bfhDQd+8QC/cuLP6nP5wwH7/PZxFOlqYvn0St1dmBEEXlCVe8WkfcouUnQv+C5ICjFmzENVKvEWJ5JOZnz03fwm5npPLN1EQd8nWFrVzZtimT8eGe7o5u/jPFeRX0er7p//x7qQIxpyMb2OY5b/5NE9gmraTrkB+K77SBzTj9ytrZg0iQrPEzdU27hoapfu3/nF68TkebA8aq6PMSxGdOgZPwQja7vR87qdiSNXc5x13zBga99bF7Yg8pf22JM7ahsh/k84AKc/+Bvgd3AfFX9VSiDqynrMPcuH+swrzr/mQsluoBmw9aROPBHOBjLr0dAv1ZWgEDdP4+1qc53mAPfuH9vxrnKCmB5ZfatCw/rMK/9fKzDvOpKm7mwSee9OvDBeeqbkKb3zPxG92bneR2m5+r6eaxN9eEmwSgRaQtcAaTVtDQzxhwr8MZCnw/+9YfmfPbAUO4cdQLvfreNsx6fz+zl24t/1BnjmcoWHn8A5gI/qOoSEekMfB+6sIxpmFJSYONGKCpy/qakQGxUJL86uzvv3TGUtk0bcdtry/jFq1+z60Cu1+GaBqxShYeqvqGq/VT1l+7yBlW9NLShGWP89WzbhLf/7zQmntOD+et2c+Zj87nnnxl07Kg24ZSpdZUqPESkm4h8LCLp7nI/d7gSY0wtioqM4BcjuvDB3cNJimzC21tXcGjwl0Q2yWHTJhg/3goQUzsq22z1IjARyAdQ5zLdK0MVlDGmfJ1aJrBxymAy5/Yhtu1+2t64gMSBP5JzSJk0yevoTENQ2Ut1l6jqKSLyjaqe5K77VlX7hzrAmrBLdb3Lxy7VDb1Ro0agKkQmHqLF2enEn7CLvG3N2DunL3PfWOZ1eCETbuexJurDpbpzgC7AMnf5MmBOZfatCw+7VLf287FLdUPP5/O/rLdI43tu0Q53zFXfvbP1yY/WaV5+odchhkS4nceaqA+X6t4GPA/0EJGtwN3ArTUq0owxNTJ5MsTHFy8JOavbs2/6CPq1aMtj/1vHBU9/xnebf/IwQhPOKnu11QZVHQ20AnoAI4GhIYzLGFOB0u4Lef6pWN6deBIvXjeQfTmHufhfi3jk/dUcOlzodbgmzJRbeIhIExGZKCJPi8hZQA4wDliPc8OgMcZDpd0XAnBWrzb871cj+PkpyTy/YAPnPLmAxRsymT7duaTXLu01NVWZUXX3AV8AtwC/BWKAi1T129CGZoypiSZx0TxySV/OP7EtE99awZUvLObQimR2b++BavSRS3vBRu01VVdRs1VnVb1eVZ8HrgIGAudZwWFM/XFal5Z8cNdwWNOJuN4ZtLtpAY067wQgJwe7tNdUS0WFR37xE1UtBH5U1azQhmSMCbZGMZFkvNuLHamnU5QXTevLl5J03jdENDpMRobX0Zn6qKJmqxNF5ID7XIBG7rIAqqpNQhqdMSZokpNh06ZmbJ82lKZD1tN0yHoaddpD5Le9UW2LiM2fbiqv3JqHqkaqahP3kaiqUX7Pa6XgEJHOIjJFRN4sb50xpnxHLu0timD/om5snzaMogPxcNo33PLKUnbst4EWTeVV9j6PahGRqSKyq3hMLL/1Y0VkrYisF5H7ykvDvUz4porWGWPKF3hpb7uERB4dfRr3/6wnn63fw1mPzWfGVxk23LuplJAWHsA0YKz/ChGJBJ4BzgF6AVeJSC8R6SsiaQGP1iGOz5gGJfDS3muvEW4e1pm5dw+nT/umTHxrBVe/+CVPTT1ol/SacoV0XktVXSAiHQNWnwqsV9UNACIyE7hQVR8BzgtlPMaY0vmSEnjtlkHMXLKZh99Zzed5C9jXujua0YlNm8Qu6TXHqNTAiDXKwCk80lS1j7t8GTBWVW92l68FBqnq7WXsnwRMBs4CXlLVR0pbV8p+44HxAG3atBkwc+bMasVvAyN6n64NjFi7rryxP0UnfU98113kbWtK5pwTyd+TSJs2ucycudjr8Ow8+vFyYMSQ1jzKUNolHWWWYKqaScA4WqWtK2W/F0RkO3B+VFTUgOoEakxDtGtjU/THgcT33E6L0Stpe/1C9n9xArsWd/E6NFOXlDViYrAeQEcg3W95CDDXb3kiMDGUMdiourWfj42qW3/5j9Yb0ShPk85bpr4Jaer7xTz9JmOf1+HZefRTH0bVDaYlQFcR6SQiMTiTSr3rQRzGmFL4j9ZbdCiGzLST2P/eQJq2KuCSfy3iT2mrbKBFE9o+DxGZgTMCb0tgJ/CQqk4RkXOBJ4BIYKqqTg5R/jYZlEf5WJ9H/fbRR6156aXO7NoVS+vWedx88wZOH7mT19ce5tPNBbRqJNzQJ5ZeSZG1HpudxxJ1fjKo+v6wZqvaz8earcLXFz/s0RF//UR9E9L00ke/U98Jh1XEae5KTQ19/nYeS3jZbBXyq628ZDUP7/Kxmkd4O1yoPPFRNCsL91N4MJa9H/bl0Po2xMYWcu+9axk9elfI8rbzWMJqHlbz8Cxfq3mY6vL5VGOO26dtb5ivvglp2vKCrzUiPld9vtDma+exhJc1Dy8u1TXGhIGMDFBtxvZ/D6XpoB9oetp64nx72PNJb1Tb2UCLYc6arSpgzVbep2vNVnXTlVcOZufOuCPL0UlZJJ2znNj2P9GvVSTjesWQ1Cj4F3TaeSxhzVbWbOVZvtZsZaorNVU1Pl6P3BMCqvEJRfp/T2zQHvfP0d4PfqCvfLFRCwuLgpqvnccSDe0+D2NMGAgcpdfngxeeF565qxMf3jOc/sc344F30rnyxcU8aQMthh1rtqqANVt5n641W9VPqsrCrQW8ml7A4QL4aWE3DizpBBpRo6uy7DyWsGYra7byLF9rtjKh5ut+SFtdvER9E9L0uOsWanSr/Qpa7auy7DyWsGYrY0zYylgXx+63B7D7nZOJanKItuM+o+mwtWRstSFO6jO7VNcYE1LO3OlCztq25G5KovmoVTQ7bT3Neu/g6039GOBr7nWIphqsz6MC1ufhfbrW51G/ffRRa/7+9+7k5ZWMg5XYbQfHX/wtORQy2hfFpV1j+Gxem2PG0yqtT8TOYwnr87A+D8/ytT4PUxtSU50+Dv8xsLJy8/XBd1aob0Ka9r//Y23afdfRl/3Glz5Wlp3HEtbnYYwJa4Fzp6ekQOPYKH5/YR/euHUI+zIjaHbRVySd8x0RsfkA5OTApEmehm3KYYWHMcZTp3RsQcaLw9j/RRcS+myl7c3zadR1B+AMgWLqJuswN8Z4Lrl9JJsW9CBnbVuSzllO60u+5uCa40hY3QeI9To8UwrrMK+AdZh7n651mIe/ozrVI4pocuoGmp3+PbFRMK5vFKe1izoy0KKdxxLWYW4d5p7lax3mpq4I7FR//KUsveRfi9Q3IU2vm/KlPj01x329qNYmnqrrbEh2Y0yDl5LiPEo05s6iIby6eBOT09bwae589rXogW7ysWkTjB9fsp+pfdZhboypsyIihHGndeTwrOHkbW1O0tkraXP1F0S1yLarsTxmhYcxps7LWB3PrtdPZc/sE4lumU27GxbSZNB6MrYUeR1ag2WFhzGmzktOBhAOpndg25Th5PzQmuYj15J84yLSt+73OrwGyQoPY0ydN3kyxMc7z4sOxrHnnQHsf/9kGrfO48JnFvG3uWuY9mqhzRlSi+xS3QrYpbrep2uX6hpwLucNHPtqyIidzFxzmIVbCyjYm8Ce9/uRt7UFQI3mDKkv7FJdu1TXs3ztUl1T35R2Hn2Ddmn7Wz/W5N+mafMz01Wi82s0Z0h9YZfqGmNMDWR81QqWDafZ8LUkDthIfNedZM7tS8bGVl6HFrasz8MYU+8lJ4PmR7Hv497snD4ELYigzRVfkXz5d/yUc9jr8MKSFR7GmHrPv0M9b2sLtr08jINLTiCi01ZGP7aAOSu2M3061qEeRNZsZYyp94rvMp80yRmJN7lDJJNv7E7/M45jwn+X88vpy8hbfxy79/RGNc7uUA8Cq3kYY8JCaXOG9G7XlHf+73T4tgcxHXfR9ub5JPTdDKjdoV5DVngYY8JaVGQEGR92YdvLw8jfnUjLc5fT+oqviGySY/OF1IAVHsaYsJecDAV7G7PztSFkftib2Hb7aHfTApLP/JHCovC91y2U6kXhISKdRWSKiLzpt+4iEXlRRGaJyNlexmeMqdtKOtSF7G86sm3KCAq2t4ABq7ji+S94fEqWdaZXUcgLDxGZKiK7RCQ9YP1YEVkrIutF5L7y0lDVDap6U8C6d1T1FuB64OdBD9wYEzZSUuCFF8DnAxHo0KIRf/3ZKTx2xYms2pLNE2s+Y1+771EpOtKZbgVI+WrjaqtpwNPAK8UrRCQSeAY4C9gCLBGRd4FI4JGA/W9U1fLGF7jfTcsYY8p07HwhAnTgV1e1IqfHSpoPX0dC9x1kzulHzs6mTJpkV2KVJ+SFh6ouEJGOAatPBdar6gYAEZkJXKiqjwDnVSZdceakfBSYo6rLghiyMaYByfg+Fl13MgdX76DF2ekcd90iDnzVmc2fd8X5PWtKUysDI7qFR5qq9nGXLwPGqurN7vK1wCBVvb2M/ZOAyTg1lZdU9RERuRMYBywBvlXV5wL2GQ+MB2jTps2AmTNnVit2GxjR+3RtYETjL9jn8corB7NzZxwAEbH5NDtjNYknbqZofyMmnaV0b1F3C5CwHxgR6Aik+y1fjlMIFC9fC/wzVPnbwIi1n48NjGhCJdjnMTVVNT5eFUoeTbvv1v73f6y+CWl6+V+Wq6/L4SNzq9eludO9HBjRq5rHEOBhVR3jLk90C7LA/o6a5mtDsnuUj9U8TKiE4jyWNtz7sJE7efzjKFYXHaAwK47MuX3J3dC6Tg313hBrHlHABqATEAN8B/QOVf5W86j9fKzmYUKlNs+jz6ca026vtr1pnvompGnSz77RiEZ5dWao97CueYjIDGAk0BLYCTykqlNE5FzgCZweqamqOjkEeVvNw6N8rOZhQqU2z+OoUSNQFYgspOmQ9TQd/ANFudHs+7gX7z2zDue6He+Efc3D64fVPGo/H6t5mFCp7ZqHf19IdKv9etx1C9U3IU1vmrZE//XyIfX51LP+kLCueXjJah7e5WM1DxMqtXkeP/qoNX//e3fy8kquuIqNy+fce5byXdE+8vMi2fdJT7KXHw9IrfeHWM3Dah6e5Ws1D1Pf1PZ5TE3VUmsXvt7Z2uaqz9U3IU1b//wLjWp6sNanvrVpaI0xpo469s50R8aqBHTlYBqfmEHzM9bQ9qb5/LSgOxnLOuHcvR7erNmqAtZs5X261mxl/NWV8+h/c2Fk4iFanJ1O/Am7KNrdhEcuymf1l8cdc/lvsJuzrNnKmq08y9earUx9U1fO47E3FxZp8/5btOfv5mrn+2ZryxHrlIjCI6/Hxwe/Q93LZqt6MSS7McbUNYEj9fp8wj/vbc/CiSMo3NiWhMHraDvuM2KO+wkg7GYutGarClizlffpWrOV8VcfzuOoUSOI67yLFmPSiUzI5cCSzuz/rBsURvDJJ/ODlo81W1mzlWf5WrOVqW/qw3ksvj9EYg5rizHL1TchTdvd8on6TtlT5tVb1WHNVsYYE0aKZy7Uw9HsnduXHTMGEREBjFrMr15bQcb2fFSp1xNPWeFhjDFBFtgfcpy05A+nD4c1nYjrnUG7mxfQqPNOoP72hVifRwWsz8P7dK3Pw/irz+dx1KgRRB+3n6RzlhPTKouDK9ux9+NeaG5MtfpCrM/D+jw8y9f6PEx9U5/P45GxsiIKtenpazX53tna4Y4P1Tdsi776alGV+0Ksz8MYYxqA4r4QiiLYv6gb26cNo+hAPJz2Lb95bylb9h6qN30hVngYY0wtCewLaZeQyKOjT4NvehLdYQ/tblpA4xMzAK3zfSE2tpUxxtSiY8fKEsZd15nIr9qQNHYFSWNXkNBzK5kf9CMjI8GrMCtkHeYVsA5z79O1DnPjLxzPY8k4WUrjfptpPmo1RBRR+G1nxg0uYOqULqWOkWUd5tZh7lm+1mFu6ptwPI+B42RFNj6kx13+lXNz4fULNbrl/lLHyLIOc2OMacAC+0I6JMXxtwsGwucnEdH4EG2v/4ymQ9dBZGGd6QuxPg9jjKkDSusLue66dsjSljQ/cyXNTv+e+G7byZzTj4yM5l6FeYTVPIwxpo5KToaiQzFkpp3EzjdOISK2gOOu/ZzkC1aRc7jA09is8DDGmDrqyH0hQO6G1mybMpzcFcnQ40fGPrGQVZmFnsVmhYcxxtRRgX0hyW2jeSylL/8ZP5jICOGvS3K577/L2X8ov9Zjs8LDGGPqsJQU2LgRioqcvykpMKhzEnPuGsa5naJ5felmzn58Pv9btbNW47L7PCpg93l4n67d52H82XkskZ2dzZ7CRkxJP8zmrCJOPS6Sa3rG0iRWgpK+3edh93nUej52n4cJFTuPJYqPxeGCQv3nx+u06+/e1xN/P1ffWrZZi4qKapw+dp+HMcaEr+jICG4f1ZX37xpK55YJ3POf77hx2hK2/XQoZHla4WGMMWHihNaJvHHraTx0fi8Wb9jLWY/NJ3XxppDkZYWHMcaEkcgI4YbTO/HhPcM5Kbk53+/MCkk+doe5McaEoeNbxPPqTaeSXxiai6Ks8DDGmDAlIsREBefKq0DWbGWMMabK6nzhISKdRWSKiLzpt66niDwnIm+KyC+9jM8YYxqikBYeIjJVRHaJSHrA+rEislZE1ovIfeWloaobVPWmgHWrVfVW4Aqg9BtYjDHGhEyoax7TgLH+K0QkEngGOAfoBVwlIr1EpK+IpAU8WpeVsIhcAHwGfBy68I0xxpQmpB3mqrpARDoGrD4VWK+qGwBEZCZwoao+ApxXhbTfBd4VkdnAa0EK2RhjTCWEfGwrt/BIU9U+7vJlwFhVvdldvhYYpKq3l7F/EjAZOAt4SVUfEZGRwCVALLBcVZ8pZb/xwHiANm3aDJg5c2a14rexrbxP18a2Mv7sPJYI6znMgY5Aut/y5TiFQPHytcA/QxmDjW1V+/nY2FYmVOw8lvByDnMv7vPYAhzvt9wB2BaKjIpH1QUOiMj31UymKbA/eFHVuXxDlU8w0w1WWi2BPUFIx3jLq89kXRTqY+Er6wUvmq2igHXAmcBWYAlwtaquDGkg1SQiL6jq+HDNN1T5BDPdYKUlIku1rCq4qTe8+kzWRV4ei1BfqjsD+ALoLiJbROQmVS0AbgfmAquB1+tqweF6L8zzDVU+wUzXq3Ng6ib7fyjh2bEI68mgjPFnNQ9jgqfO32FuTBC94HUAxoQLq3kYY4ypMqt5GGOMqTIrPIwxxlSZFR5hRkQSROTfIvKiiKR4HY8xdUVpI3Q3VME4FlZ4hIiINHOHjF8jIqtFZEg10yl1ZGL3tdJGJ74EeFNVbwEuqMFbCGv2ReINETleRD51PxMrReSuGqRVpc+GljJCt5dEJE5EvhKR79xj8fsapFXrx8IKj9B5EvhAVXsAJ+Lc03KEiLQWkcSAdSeUks40AkYmdrctdXRinDv2N7ubFdbwPdQrVZkCoK59kTQgBcCvVbUnMBi4zf2/PSKEn426Jg8YpaonAv2BsSIy2H+DunwsrPAIARFpAgwHpgCo6mFV/SlgsxHALBGJc/e5BXgqMC1VXQDsLSWbI6MTq+phYCZwIc7wLx3cbRra+Z1GJacAqP3QDICqblfVZe7zLJwfVe0DNgvVZ6NOcYePynYXo91H4OWvdfZYNLQvl9rSGdgNvCwi34jISyKS4L+Bqr4BfADMdPsmbsSZ3Kqy2lNSwwCn0GgPvAVcKiLP0sDuxC3jA1QvvkgaInfoopOAL/3Xh+qzISJJIvIccJKITKxJ7MEiIpEi8i2wC/ifqtabY2GFR2hEAScDz6rqScBB4JgZE1X1r0Au8Cxwgd+vkMoobVZ7VdWDqnqDqv5SVadXI/ZwU2++SBoSEWkM/Be4W1UPBL4eos9Gpqreqqpd1Jk/yHOqWqiq/XFaC04VkT6lbFMnj4UVHqGxBdji9yviTZzC5CgiMgzoA7wNPFSNPGpldOJ6rt58kTQUIhKNU3BMV9W3ytimQX023GbteZTeb1Enj4UVHiGgqjuAzSLS3V11JrDKfxsROQl4EacJ5QaghYj8qQrZLAG6ikgnEYkBrgTerXHw4afefZGEMxERnL7A1ar6WBnbNIjPhoi0EpFm7vNGwGhgTcA2dfdYlDXRhz1qPAlWf2ApsBx4B2ge8PrpQF+/5WjgllLSmQFsB/Jxvghv8nvtXJzh7X8AJnn9nuvCg2MnH4sCNgCdgBjgO6C313E21AcwFKdTeDnwrfs4N2CbBvHZAPoB37jHIh14sJRt6uyxsLGtTNgQZwqAkTiTPu0EHlLVKSJyLvAEEAlMVdXJngVpTJiwwsMYY0yVWZ+HMcaYKrPCwxhjTJVZ4WGMMabKrPAwxhhTZVZ4GGOMqTIrPIwxxlSZFR4mbImIisg//JbvFZGHg5R2VcYXqjNpVzL/ju6xu8Nv3dMicr2HYZk6xgoPE87ygEtEpKUXmYtIlBf5Bsku4C53SAtjjmGFhwlnBcALwD2BL4iIT0Q+FpHl7t9kd/00EXlWnNnuNojICHeSqdUiMi0gjX+IyDJ3/1buunki8mcRmY/z5TtAROaLyNciMldE2pYSSycR+UJElojIHwNe+427frmUMdOciGSXEcst7r7fich/RSTeXX+5iKS76xeUcex2Ax8D48o9wqbBssLDhLtngBQRaRqw/mngFVXtB0zn6Al2mgOjcAqd94DHgd5AXxHp726TACxT1ZOB+Rw92mkzVR3hpvlP4DJVHQBMBUobGuVJnOH7TwF2FK8UkbOBrjhzkvQHBojI8FL2LyuWt1T1FHVmqlsNFM+c+CAwxl1f3lTFjwK/difUMuYoVniYsKbOXBGvAHcGvDQEeM19/irOgH3F3lNn3J4VwE5VXaGqRcBKnIEXAYqA/7jPUwP2L17fHWco7f+5E/7cT8ksj/5OxxnYrjiWYme7j2+AZUAPnMIkUFmx9BGRhSKyAkjBKQABFgHT3FnpyiwYVPVH4Cvg6rK2MQ1XfW6TNaaynsD58n25nG38B3nLc/8W+T0vXi7rM+O//0H3rwArVXVIJWIsbZA5AR5R1ecrsX9paU0DLlLV79zO7pEAqnqriAwCfgZ8KyL9VTWzjLT+jDMfTVnNW6aBspqHCXuquhd4nZJmG4DPceY2AOdX+WdVTDYCuMx9fnUZ+68FWonIEHAmQRKR3qVstygglmJzgRvdWfcQkfYi0roKsSQC293Jl46kKyJdVPVLVX0Q2MPR850cRVXX4MxFc15Z25iGyWoepqH4B3C73/KdwFQR+Q1O5/ANVUzvINBbRL4G9gM/D9xAVQ+LyGXAU26fSxROLWhlwKZ3Aa+JyF04M+wV7/+hiPQEvnDmUCIbuAbnSqjKxPIAzvzgm3Ca4BLd9X8Tka44NZuPceY4Kc9knKYzY46wIdmNqedEJFtVG3sdh2lYrNnKGGNMlVnNwxhjTJVZzcMYY0yVWeFhjDGmyqzwMMYYU2VWeBhjjKkyKzyMMcZUmRUexhhjquz/Ac3edzAbddZ+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 4.\n",
    "\n",
    "plt.loglog(N_array, residuum, 'bo', label='Residuum')\n",
    "plt.title('Convergence de la méthode de Simpson')\n",
    "plt.grid(which = 'both')\n",
    "plt.ylabel('Residus')\n",
    "plt.xlabel('Nombre de pas N')\n",
    "\n",
    "TN = np.logspace(0.7, 1.5)\n",
    "plt.loglog(TN, residuum[0]*(N_array[0]/TN)**6, label='$1/N^6$')\n",
    "\n",
    "plt.legend()\n",
    "#gradient = np.polyfit(np.log(N_array), np.log(residuum), 1)[0]\n",
    "#print(f\"L'intégrale converge en res(N)~ N^({gradient:0.1f})\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}