Cuadernos/MN_Butcher.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# La Tabla del Carnicero (Butcher's Table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introducción"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ya hemos introducido dos métodos numéricos para resolver ecuaciones diferenciales, a saber el método de **Euler** y el de **Runge-Kutta*.\n",
    "\n",
    "Se preguntará el espectador, ¿Existen más?\n",
    "\n",
    "**¡Por suspuesto que sí!**\n",
    "\n",
    "Hoy veremos una forma de construir muchos más mediante un proceso mnemotécino llamado la **Tabla de Butcher**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**No olvidemos nuestro foco de atención**,\n",
    "\n",
    "$$\\frac{dx}{dt}=f(t,x),\\text{     }x(t_{0})=x_{0}$$\n",
    "\n",
    "esto es, la tasa de cambio instantáneo de la cantidad dependiente $x$, se puede modelar mediante la función $f(t,x)$. Además el experimento comienza en el momento $t_0$ con una cantidad inicial $x_0$ de la cantidad $x$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definición"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La tabla del carnicero es una regla mnemotécnica que ayuda a construir métodos numéricos ya conocidos mediante una tabla que resume su información. \n",
    "\n",
    "El objetivo de esta lección será acostumbrarnos a la notación de tabla para construir y programar el respectivo método, así como el proceso inverso, esto es, ir del método numérico a la tabla.\n",
    "\n",
    "Esto es muy importante, pues el entender la tabla genera un proceso de abstracción importante en nuestro cerebro que nos hará entender la estrucutra interna de los llamados **métodos numéricos explícitos** para resolver ecuaciones diferenciales.\n",
    "\n",
    "El proceso se puede resumir con las siguiente gráficas (Tomadas de [Wikipedia](https://wikimedia.org/api/rest_v1/media/math/render/svg/44d56b38e104bfdf185d1307337e5c6e074b16c9)):\n",
    "\n",
    "**Ir de**\n",
    "\n",
    "![Tabla](/files/Pictures/Butcher.png)\n",
    "\n",
    "\n",
    "**Para construir el método**\n",
    "\n",
    "![MetNum](https://wikimedia.org/api/rest_v1/media/math/render/svg/44d56b38e104bfdf185d1307337e5c6e074b16c9)\n",
    "\n",
    "![Predictor](https://wikimedia.org/api/rest_v1/media/math/render/svg/84fa47cf3289c5ab2067fe46f748f33d4b6be48c)\n",
    "\n",
    "**y vicerversa. Recuerde que en nuestro problema, $x$ hace el papel de la $y$ mostrada arriba.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aunque esto puede ser un poco confuso al principio, iremos paso a paso a paso para ver cómo se hace dicha transformación.\n",
    "\n",
    "- **c's:** Representan dos cosas:\n",
    "\n",
    "    1. Cuántos **k's** necesitamos definir para construir todo el método numérico.\n",
    "    2. Cuánto se le va a sumar a $t_i$ en multiplcación con $h$ en $f(t_i,\\dots)$ **antes de la coma**.\n",
    "    \n",
    "   \n",
    "- **a's:** Representan cuánto de cada k se le debe agregar a la función en multiplicación con $h$ $f(\\dots,x_i)$ **después de la coma**.\n",
    "\n",
    "- **b's:** Ayudan a construir la función **predictora** del método (futuro = presente + paso\\*( comb. lineal de **b's** y **k's**). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ejemplos básicos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Método de Euler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Tabla](/files/Pictures/BEuler.png)\n",
    "\n",
    "Como podemos ver, en este método, necesitamos una k y sumar $hk_1$ en el predictor.\n",
    "\n",
    "El método se escribe entonces, como:\n",
    "\n",
    "$$k_1=f(t_i,x_i)$$\n",
    "$$x_{i+1}=x_i + hk_1$$\n",
    "\n",
    "que si recordamos, es lo mismo que tener:\n",
    "\n",
    "$$x_{i+1}=x_i + hf(t_i,x_i)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Método del Punto Medio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Tabla](/files/Pictures/BMidPt.png)\n",
    "\n",
    "En este método, necesitamos dos k's y sumar $h\\bigl(\\frac{1}{2}k_1+\\frac{1}{2}k_2\\bigr)$ en el predictor.\n",
    "\n",
    "El método se escribe entonces, como:\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "k_1 &=f(t_i,x_i) \\\\\n",
    "k_2 &=f\\biggl(t_i+\\frac{h}{2},x_i+h\\frac{k_1}{2}\\biggr) \\\\\n",
    "x_{i+1} &=x_i + h\\frac{k_2}{2} \\\\\n",
    "\\end{split}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Método de Ralston"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Tabla](/files/Pictures/BRalston.png)\n",
    "\n",
    "En este método, necesitamos dos k's y sumar $h\\bigl(\\frac{1}{4}k_1+\\frac{3}{4}k_2\\bigr)$ en el predictor.\n",
    "\n",
    "El método se escribe entonces, como:\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "k_1 &=f(t_i,x_i) \\\\\n",
    "k_2 &=f\\biggl(t_i+\\frac{2h}{3},x_i+h\\frac{2k_1}{3}\\biggr) \\\\\n",
    "x_{i+1} &=x_i + h\\biggl(\\frac{1}{4}k_1+\\frac{3}{4}k_2\\biggr) \\\\\n",
    "\\end{split}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora veremos un caso curioso, llamado **Runge-Kutta 3/8**:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementación"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Haremos la implementación de Runge-Kutta, Regla 3/8.\n",
    "\n",
    "Los demás métodos mostrados y existentes son deber el estudiante.\n",
    "\n",
    "**Por favor, documentar adecuadamente el siguiente código:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RK38(f,x0,t0,tf,n):\n",
    "    t=np.linspace(t0,tf,n)\n",
    "    h=abs(t[1]-t[0])\n",
    "    x=np.empty(n)\n",
    "    x[0]=x0\n",
    "    \n",
    "    for i in range(n-1):\n",
    "        k1=f(t[i],x[i])\n",
    "        k2=f(t[i]+h/3  ,x[i]+h*k1/3)\n",
    "        k3=f(t[i]+2*h/3,x[i]+h*(-k1/3+k2))\n",
    "        k4=f(t[i]+h    ,x[i]+h*(k1-k2+k3))\n",
    "        \n",
    "        x[i+1]=x[i]+h*(k1/8+3*k2/8+3*k3/8+k4/8)\n",
    "    return t,x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Es claro que esta debe ser la función de un módulo, junto con ```EulerM```, ```RK4``` y todos los demás métodos.\n",
    "\n",
    "**Observación:** El método RK38 contiene meno error que el método clásico RK4 (aunque requiere un poquito más de cálculos, pero la diferencia en tiempo de cómputo es insignificante).\n",
    "\n",
    "Resolvamos la siguiente ecuación con la regla 3/8:\n",
    "\n",
    "$$\\frac{dx}{dt}=\\tan(x)+1,\\text{    }x(1)=1$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtwAAAGLCAYAAAD0ykW7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3hVVdbH8e9KAgmkQRI6ISBIMfQqZRQsCCooFhQBUextZJDXGRswgjo4lnEKMhbEgljGShEVNShSBAuMqFSRQJQqkIQkkGS/f5ybTAipmOSm/D7Pcx9yzz5lnbMDrLvvOvuYcw4RERERESkfAf4OQERERESkOlPCLSIiIiJSjpRwi4iIiIiUIyXcIiIiIiLlSAm3iIiIiEg5UsItIiIiIlKOlHCLSKVlZo+Y2TdmVsffsYiIiJwoJdwiUi7MbKqZOTPbVEj7Zl/71ELaBwDjgIucc2l5lvcubJvfEOsaM5tTlvvMs++OvvMcWAb72ubblzOzI2a2ycxmmFloGYTqN77zuTXfst5mdtDMFptZcCn2dWf+a21mtX2/j13LKOSc/V7s+z0O9L0fmKd/nJkdMLNVZnZhAdtuM7NH8i07ycx2mNlXZlavgG3eLuRaLTSz+8ry3ESkbCnhFpHylA60MrOeeReaWS8gztd+HF8CORsY45zbmq+5NzClHGKtKl4G+gJnAS8AfwCe8GtEZcyXGC8GvgRGOOcySrH5ncDAfMtq4/3OlFnCbWYBwJ+BvzrnsvI1j8broyuAfcCbZnZaMftrAXwM/AoMds4dyNc+GDi1kM3/AkwsKEkXkcpBCbeIlKdUvCTi8nzLL/ctTy1oI+dcqnOurXPu/XKOryr62Tm30jn3qXNuGr4PJr4EsMozs3jgQ+A7YFjebzcqmTOB1ngfgPJb5+ujRcAleEn0mMJ2ZGZNgY/wPoCe5Zzbm6+9Ft6HqnsK2t459xleYj/2BM5DRCpAtfgHWkQqtVeAkWZmAL4/R/qWH8fMBpjZUjM7bGb7zOxpMwv3tV0F/MP3c87X9gl5tj3D9xV+upntMrOZZhaWb/8dzexz3zrfm9nwQuIYaWb/NbMMM0s0swfMLKi4kzWzm33rp5rZfKBJAesEmNmffOUIGWa20czGFbfvQqwFgoEGefafYGb/yXfMnHKHjr73LX3vR5rZv33lGzvM7M/5k3czu9RXvpJmZp+YWTfftlflW+9aM1vvO6efzOzO0pyImbUFlgBbgXOdc6l52qaa2d4CtsktsTCzbUA0MCXP78dAINm3+nN5lrf0bfMXXz+n+M5/rpk1LkG444APnHPJRa3knDsMbAZiCznnhnjJtgFnOOd2FbDa7UAa8FwRh3oDuLIEcYuIHyjhFpHy9ibQCBjge/87vOTwrfwrmll/vOTjF7yRwQnAufwv0VgIPOr7ua/vdbNv21PwyhD2AhfjlRBcAeQmnubdfPk+EOZrmw78DWiRL47BwKvAV8AFeEn+JOCfRZ2omV0A/AtYAFwE/BdvBDq/fwD3Ak8B5/muxWwzO7+o/ReiBV5CeVwyWkIPAyl41/slYLLvZwB85UCv4F2LEcC7eNfmGGb2f8CTwNvA+b6fp+WvNy5CK7y+/xk4xzl36ATOZQRwEHiW//1+fAWc4Wufnmf5z75lDYEH8fphAnAS8LH56rKLcAawvLiAfB9emgM/FtAchfcBow5esp1UwPaNgfuACc657CIOtRzoYWb1i4tJRCpesaM1IiK/hXPugJktxisj+cz352Lf8vyr/wVY7py7LGeBme0EPjKzjs65b32jmDjnVubbdjLwEzA8p6bWzPYDr5pZX+fcCuBqvASrj3Nuh2+dbcCyfPu6H0hwzuWMOi/2xfqQmU3P2bYA9/jO7Sbf+/fNrAFwbZ7zaQPcBFztnHvet3iJmTXB+5CwoJB959mFBeHVJZ8G3Ag8UEAdcUl96py7w/fzh2Y2BO/Dwmu+ZX8Evgcud845vGtRC5iRJ6AIX+zTnXN/zrOvusC9ZvZkCeKbiDeK2yt//XJJOee+NrNMYEfe3w8zW+37cUv+3xvn3Pg86wUCK4AdQH/g04KO4ysBaQJ8W0gogb4+qo9XUx5KwXX2VwMO6Oac217Ivh4G3nfOFRhLHmvxRsl74pXkiEglohFuEakIrwCXmDfbxCUUUE7iS876Aq+ZWVDOCy8ZPgr0KOYYvYG38iV2bwCZ/G90vTfwZd6E2Tn3ObA7TxyBQHfg9Xz7fxXv38y+BR3ct1034J18TW/me38mkA28le88PwK6lmBkdSLe9UgF3gM+cc7NKHqTIn2Q7/13eCOyOXoB833Jdo53823TFy+pfD3fOX2M9+1Gc/ME5XnlP88leKUxf7ECPomVFzMbambLzewg3u9Kzu9G2yI2yyk5KexbhW/w+mg3Xn9d5ZzbUMB6n+HVbT9iZrULiK0v3t+X/yv2RP4XS0nKYUSkginhFpGK8C5eGccDeInZ/ALWqQ8EAjPxkpWcVwZQi0JqYPNoAhxT/+pLvvfhfXUPXjKym+PlXRbjO17+Wtqc91EUrAHet4b595//fQzeeR7k2POc49v+uJrvfF7CS4IH4pXajDCzm4rcomj5R5OPACF53jcG9uRbJ//7GN+f6zn2nD7xLY8FTs/X9lG+fbwD3IZXG/1bPkCUmHmz5byLl2SPxfvgkDMTSEhh2+VpK2z2lMvx+uhiYANe7XjTAtb7wrfuIOCFAj5o/A34N3DQzOrZ/2YhqWNmkfnWzYmlqLhFxE9UUiIi5c45l2pmC/CmsHs9781weRzA+3p9KrCogPbj6lvz+RmvXCSXbxQ1GtjvW/QL0L6AbfNutxcvIWyYb51Gvj/3U7A9eCOk+bfL/36/b73+eCPd+RX0gSCvXc65Nb6fl5pZHHC/mb2Q57qm45Wc5FXYB4Xi/EKeGzJ98r/PuSbnc/wHFfCSTvCS0BzH3WzonJuZU7NsZr845x7L03zcOZVBvfIIvH67LGcE33c9i5NzvoVNw7feOfctsMbM1uJ9a3AfXinRMZxz75rZ9Xh157uB3+dpbof3rcyEfJs9DDzEsf+H58RS2O+niPiREm4RqShP4pUMzCqo0ZeUrwTaOefuL2I/RwDMLMQ5l3ce71V4o7135ykruQjv37mcGu3VwGgza56nhrs/eZJi51yWmX0JXOqLOcdIvAR5RSHxZ5nZN3g3WeY9x4vyrfox3gh3pHOuLGpt78I792uAv/uW7cCr787r7BPc/2pgmO+65pSV5J/ZZQVe/XVT59zCIva1pog2AJxzk30zdzxiZrudcy/5mnYA4WbWzDm307dscAG7yD9Cn7OMApbXAY7mK5cZXVyMeDdAHsG70TOhqBWdc1vM7BlgvJlNcc4d94HKOTfbzBoBD5rZLufcA76m8zn+/+lP8Po5f6lSS9+fG0sQv4hUMCXcIlIhnHMJFJOc4N1g9pGZZePNLpKMNwvHecA9zrmNwA++dW83s4+BQ7762OnA18DbZvYkXh3yDLwbznKS5OfwZgdZaN7TKusA0zi+FncK3g2Pz+HVm3fyrfd0ETdMgjfbxZu+47+FV0YxJN912GBms4BXzOxhvCQ0BIgH2jrnrqUUnHNfmNmHwB/M7F++DxtvAdeY2eN4M7sMAs4pzX7zmIGX0L/iux4dgOt8bdm+GA74rucTvhHiT/FKFtsCg5xzI0p5zJvxPgTNNrO9zrnFeDPQpPmWPYqX7N5YwLY/AOf5btRNATY455LN7Ee86Sm/xRstX4d3c+EEM/sbXplTP4qYLzuHcy7D96GsB0VP1ZfjYbxrdhveSHdB+3zIN7o/3Te6/6xzLv/NvPiqTjY555bma+qJV6a0vgTxiEhFc87ppZdeepX5C680ZG8x6+wFpuZb1gcvuTqEd2Pgd8BjeCPC4M3E8DBeiUk23mwiOdueiZccpuN9PT8TCMu3/854U6hl4JU6XIiX9M7Jt95leNP6HcEbXX0ACCrBed/qW/8wXmnMYLxSmYF51jG8MoH1vjj2AEuBK4vZ9zbgkQKWn+Y7xqg8y+4CEvE+tLyENyrtgI6+9pa+9+fn29ccYE2+ZSPx5pJOx/u24CzfthfmW28M3tMh0/Ae9rIKmFiCa+aAW/MtC8b7gJaCN6sMwFDfNTuMd8Nhh/zb4iXBK32/O7nX3dcP63zn4ICWvuV3+q5TKt6NmycXFE8BMU8CNudbNjDvNc7XNhvvfoLQwvrS93vxMl7J0QUlvVa+5e8Az/n7771eeulV8Mucy/tNmoiISNHMbAzwInCSc66g+aWrPV8JyHZggHNudXHrl3MskXi182e5AkbFRcT/lHCLiEiRfCUyH+KNWnfHK8v53Dl3Ig/qqTbM7F9437wUW4ZSznH8CRjinBvozzhEpHCq4RYRkeJE45XnROOVRbyKV4pR003Dq5UPdCf+4KGycJBjZzcRkUpGI9wiIiIiIuVID74RERERESlHSrhFRERERMpRhddwm9lsvMn8dzvnOhbQfgFeXVw23tRIE3LuujazLLxpugC2O+fyP3zhODExMa5ly5ZlFH3ppKamEhoa6pdjS8VQH9cM6ueaQf1c/amPawZ/9vOXX3651zmX/2m8FV/DbWan4c2r+kIhCXcYkOqcc2bWGXjNOdfe15binAsrzfF69uzp1qwp9uFm5SIhIYGBAwf65dhSMdTHNYP6uWZQP1d/6uOawZ/9bGZfOud65l9e4SUlzrlPgf1FtKe4/30KCMWb5F9EREREpEqqlDXcZjbCzH7AeyTx+DxNIWa2xsxWmtmFfgpPRERERKTE/DItoJm1BBYUVFKSb73TgMnOubN875s655LM7CTgY+BM59yWAra7HrgeoFGjRj1eeeWVMj6DkklJSSEsrFQVMFLFqI9rBvVzzaB+rv7UxzWDP/t50KBBBZaUVOoH3zjnPjWz1mYW45zb65xL8i3famYJQDfguITbOfcU8BR4Ndz+quNRrVj1pz6uGdTPNYP6ufpTH9cMlbGfK13CbWZtgC2+mya7A7WBfWZWHzjsnMswsxigP/DwiR7n6NGj7Nixg/T09LIJvACRkZF8//335bZ/8T/1cc2Q08+BgYHUq1ePmJgYAgIqZUWeiIhUQv6YFnAeMBCIMbMdwBSgFoBzbhZwMXClmR0F0oDLfMl3B+DfZpaNV3v+F+fcdycax44dOwgPD6dly5aY2W87qUIkJycTHh5eLvuWykF9XDMkJycTFhbG0aNH2bVrFzt27KBFixb+DktERKqICk+4nXOjimmfAcwoYPlyoFNZxZGenl6uybaIVC9mRu3atWnWrBkbNmzwdzgiIlKF1OjvRJVsi0hpqZRERERKS/9ziIiIiIiUIyXc1ci2bdswMzIzM3/Tfm688UamTZt23PLXXnuNc845h4yMjN+0f3/4/PPP6dWrF/v3F/rMJQDi4+NJSEiomKBERESkRlDCXQktW7aMfv36ERkZSVRUFP3792f16tUVdvxZs2Zx3333HbPs66+/5tlnn+Wtt94iODj4hPY7depUzIzXX389d1lmZiZmxrZt235LyEVKTEzk7rvvZuHChURFRRW57vr16yvdVEIlMWfOHAIDAwkLCyMiIoIuXbqwYMGC3Pb8H8acc9x22220b9+enTt38sknn9CpUyfq1atHdHQ0I0aMYOfOnbnb79+/n8suu4yYmBhiYmIYPXo0hw4dqvDzLMp9991Hp06dCAoKYurUqce1v/zyy8TFxREaGsqFF154zIev/fv3M2LECEJDQ4mLi+Pll1+uwMhFRKS6U8JdyRw6dIjzzz+f2267jf3797Nz506mTJlywkluWenWrRvvv/8+devW/U37iYqKYvLkyWRlZZVRZMWLjY1l6dKlNGzYsNB1fuu3ApVB3759SUlJ4cCBA9x8881cfvnlHDhw4Lj1nHPccMMNJCQksHTpUpo1a8Ypp5zC+++/z4EDB0hKSuLkk0/mpptuyt3m3nvv5ddff2Xr1q1s2bKFXbt2FZjU+lObNm14+OGHOe+8845rW79+PTfccAMvvvgiu3btom7dutx888257bfccgu1a9dm165dzJ07l5tuuon169dXZPgiIlIGjmZl+zuEAinhrmQ2btwIwKhRowgMDKROnToMHjyYzp07A5Cdnc306dOJi4ujYcOGXHnllRw8eLDAfbVs2ZIlS5bkvp86dSpjxozJfZ8zkl6vXj1iY2OZM2cOAFdddRX33ntv7npPP/00bdq0ISoqiuHDh5OUlJTbZmbMmjWLk08+mfr163PLLbdQ1NNLhwwZQu3atXnppZcKbB84cCDPPPNM7vs5c+YwYMCAY443c+ZMTj75ZMLDw7nvvvvYsmULffv2JSIigpEjR3LkyJHc9RcsWEDXrl2pV68e/fr1Y926dcdcnxkzZtC5c2dCQ0PJzMw85pplZWXx4IMP0rp1a8LDw+nRoweJiYkA3H777cTGxtKsWTN69OjBZ599Vug5p6WlcccddxAXF0dkZCQDBgwgLS0NgHfffZf4+Hjq1avHwIEDj5nTu2XLljzyyCN07tyZyMhILrvsshLNGx8QEMDYsWNJTU1l06ZNx7RlZWVx1VVXsWbNGhISEmjUqBEAjRo1omnTprnrBQYGsnnz5tz3P/74IxdeeCERERFERkYyYsSIEiekOX04adIk6tevT6tWrXjvvfdKtG1pjBs3jqFDhxY4TePcuXMZNmwYp512GmFhYUybNo0333yT5ORkUlNTeeONN5g2bRphYWEMGDCA4cOH8+KLL5Z5jCIiUn7WJh7grMeW8u3eyjeIVukefONPl/17xXHLzu/chLF9W5J2JIurnvviuPZLejTn0p6x7E89wk0vfXlMW1ZWFuP6n8SwLk1JOpBG03p1io2hbdu2BAYGMm7cOC6//HJOPfVU6tevn9s+Z84c5syZwyeffJKbcN96662lTg62b9/O0KFDeeqpp7jkkks4dOhQbjKZ18cff8xdd93FBx98QHx8PJMmTeLyyy/n008/zV1nwYIFrF69mkOHDtGjRw+GDRvGkCFDCjyumTFt2jQmTJjAFVdccUIzxSxevJgvv/ySxMREunfvzvLly5k7dy7R0dH07duXefPmMW7cOL766ivGjx/P/Pnz6dmzJy+99BLDhw9nw4YNud8YzJs3j4ULFxITE0NQ0LF/HR577DHmzZvHokWLaNu2LevWrcsd4e/VqxeTJ08mICCA2bNnc+mll7Jt2zZCQkKOi3fSpEmsX7+e5cuX07hxY1atWkVAQAAbN25k1KhRvP322wwcOJDHH3+cYcOG8d1331G7dm3Aq5tfvHgxISEh9O/fnzlz5nDjjTcWeX2ysrJ47rnnqFWrFnFxcce0jR49msTERD7++GPq1at3TNv27dvp3Lkzhw4dIjAwkKeffjq37ZZbbmHmzJmMGuXN6vnGG28wfPjwknQXAKtWrWLcuHHs3buXp556imuuuYadO3cW2P/nn38+y5YtK3A/AwYMOKZUpqTWr19Pv379ct+3bt2a2rVrs3HjRgICAggMDKRt27a57V26dGHp0qWlPo6IiPjHnM9/5IFF39MgLJiQwMo3C51GuCuZiIgIli1bhplx3XXX0aBBA4YPH86uXbsAb6Ru4sSJnHTSSYSFhfHQQw/xyiuvlLokYu7cuZx11lmMGjWKWrVqER0dTdeuXQtcb/z48XTv3p3g4GAeeughVqxYcUzN9Z/+9Cfq1atHixYtGDRoEN98802Rxx4+fDgNGjQ4ZiS7NP74xz8SERFBfHw8HTt2ZPDgwZx00klERkYydOhQvv76a8Abmb/hhhvo06dP7oeY4OBgVq5cmbuv3//+98TGxlKnzvEfhp555hmmT59Ou3btMDO6dOlCdHQ0AGPGjCE6OpqgoCDuuOMOMjIyCpybOTs7m9mzZ/PEE0/QrFkzAgMD6devH8HBwbz66qucd955nH322dSqVYtJkyaRlpbG8uXLj4mvadOmREVFMWzYsCKv7cqVK6lXrx4hISFMmjSJl1566bgymg8++ICRI0cel2wDtGjRggMHDrB3716mT59O+/btc9u6d+/OkSNHiI6OJjo6msDAwGNKMooTFxfHddddl9sPP//8c+7vdH4LFizgwIEDBb5OJNkGSElJITIy8phlkZGRJCcnF9kmIiJVw67kDE5v24BFt/+ONvUD/R3OcTTCncerN/QttK1O7cAi26NCax/XnvcphCUZ3c7RoUOH3PKOH374gTFjxjBhwgTmzZtHUlLSMaOWcXFxZGZmFpq8FCYxMZHWrVsXu15SUhLdu3fPfR8WFkZ0dDQ7d+6kZcuWADRu3Di3vW7duqSkpBS73+nTp3P11VczduzYUsUN5JZBANSpU+e497/88gsAP/30E88//zz/+Mc/ctuPHDlyTElMbGxsoccp6ho9+uijPPPMMyQlJWFmHDp0iL179x633t69e0lPTy9wP/n7MiAggNjY2GNuVsx/bfPGnt+pp57KsmXLSElJ4ZprruGzzz5j5MiRx6yzYMECzj//fOrXr8/48eML3E9UVBTjxo2jS5cu7Ny5k6CgIC699FK6dOnCO++8g3OOSZMmMWbMGF577bVC48kr/3kAJfo9KSthYWHH3eR56NAhwsPDCQgIKLRNREQqr6+2/0p2tqNnyygmDW5HgFXeZ6xohLuSa9++PVdddRXffvstAE2bNuWnn37Kbd++fTtBQUHHJJ05QkNDOXz4cO77nEQUvERzy5YtxR4///FSU1PZt28fzZo1O6HzyXH22WfTpk0bZs6cWeKYSys2NpZ77rnnmBHSw4cP55ZFQNF/MQu7Rp999hkzZszgtddeY/v27Rw4cIDIyMgCa9djYmIICQkpcD/5r61zjsTExN98bcPCwpg5cyYvvvhi7mh/jn79+jF//nxuv/32ImfiyMzMZPfu3bmJ6Nq1a7nhhhsIDQ0lLCyMG2+8kUWLFv2mOAszdOhQwsLCCnwNHTr0hPYZHx/P2rVrc99v3bqVjIwM2rZtS9u2bcnMzDym3n3t2rXEx8f/5nMREZGyl53t+PfSLYyctYIZi3/AOUdggFXaZBuUcFc6P/zwA48++ig7duwAvFHWefPmceqppwLezZSPP/44P/74IykpKdx9991cdtllx9UfA3Tt2pVXXnmFo0ePsmbNGv7zn//kto0ePZolS5bw2muvkZmZyb59+wosV7jiiit47rnn+Oabb8jIyODuu++mT58+uaPbv8UDDzzAww8/fFzMb775JocPH2bz5s08++yzJ7z/6667jlmzZrFq1Sqcc6SmprJw4cISlwpce+213HfffWzatAnnHOvWrWPfvn0kJycTFBREgwYNyMzM5P777y90iryAgADGjx/PxIkTSUpKIisrixUrVpCRkcHIkSNZuHAhH330EUePHuXRRx8lODj4mFrjExUdHc21117L/ffff1zb6aefzptvvsn111+f+zvx5ptvsmHDBrKzs9mzZw8TJ06kW7duudMo9urVi2eeeYa0tDTS0tJ46qmn6NKlS+4+Bw4cWGazlrz33nukpKQU+CrqZsujR4+Snp5OdnY2mZmZpKen586GM3r0aObPn89nn31GamoqkydP5qKLLiI8PJzQ0FAuuugiJk+eTGpqKp9//jnvvPPOCX37IiIi5WtPcgZXzVnNQ+/9wNmnNOKZcb0qdaKdQwl3JRMeHs6qVavo06cPoaGhnHrqqXTs2JFHH30UgPHjxzN27FhOO+00WrVqRUhIyDElE3lNmzaNLVu2UL9+faZMmcIVV1yR29aiRQsWLVrEo48+SlRUFF27dj1mBDDHmWeeybRp07j44otp0qQJW7Zs4ZVXXimTc+3fvz+9e/c+Ztkf/vAHateuTaNGjRg3bhyjR48+4f337NmTp59+mltvvZX69evTpk2b3FKdkpg4cSIjR45k8ODBREREcM0115CWlsY555zD0KFDadu2LfHx8YSEhBRZmvLII4/QqVMnevXqRVRUFH/84x/Jzs6mXbt2vPTSS9x2223ExMQwf/585s+fn3vD5G81YcIEFi1adMzMLDnOPvtsXn31Va666irmz5/Pzp07GTJkCOHh4XTq1ImAgADeeuut3PVnz57Ntm3baN68Oc2aNWPr1q3HXMvExET69+9fJnGfqOuuu446deowb948HnjgAerUqZN7M3F8fDyzZs1i9OjRNGzYkOTk5GO+XZk5cyZpaWk0bNiQUaNG8eSTT2qEW0Skkkncf5ihT3zGqq37mH5hR2aO7k5knVr+DqtErKgp3KqDnj17ujVr1hy3/Pvvv6dDhw7leuy8NdxSPamPYceOHVx66aWsWHH8LD/VRf5+roh/P6TiJSQkVMkHX0nJqY+rtuxsx/0LvmNU7xa0a1z4/73+7Gcz+9I51zP/co1wi8hv0rx582qdbIuIiP8k7j/MVc99wc4DaQQEGFOHxxeZbFdWmqVERERERCqdd77Zyb1veZNG/LgnlWalmPGtslHCLSIiIiKVRkpGJpPf+ZY3v9pJ9xb1eOLybsRG1fV3WL+JEm4RERERqTSeWLKRt7/eye1nnsxtZ7QhKLDqV0Ar4RYRERERv8rKduxPPUKD8GB+f+bJDOnYmB5xUf4Oq8xU/Y8MIiIiIlJlJR1IY/QzKxn77CqOZGYTHlKrWiXboBFuEREREfGT+WuTuOet/5KV7ZgyPJ5agZX/ITYnQgm3iIiIiFSo1IxM7nv7W978eifdWtTjb5d1JS461N9hlRuVlNRwc+bMYcCAAbnvw8LC2Lp1a5HbbN++nbCwsNzHZudIS0ujf//+LF68uFxizS8+Pp6EhIRi1yvJOQFs27YNMyMzM7MMohMREZHCBAUam/ekcPuZJ/P6DX2rdbINSrgrrZdffpmePXsSFhZGkyZNGDp0KMuWLSv346akpHDSSScVuU6LFi1ISUkhMDDwmOU33HADkyZNYsiQISd8fDNj8+bNJVp3/fr1JXqSVEnOSURERMrX0axsZi3dwsG0owQHBfLGTf34w9ltq8UsJMWp/mdYHubOhZYtISDA+3Pu3DLd/WOPPcaECRO4++672bVrF9u3b+fmm2/mnXfeKdPjlLUXXniBESNG+DsMERERqaFjVmgAACAASURBVGQ2707hopnL+ct7P7Bw3c8A1KoBiXaOmnOmZWXuXLj+evjpJ3DO+/P668ss6T548CCTJ0/mX//6FxdddBGhoaHUqlWLYcOG8de//hWAjIwMJkyYQNOmTWnatCkTJkwgIyMDgISEBJo3b86jjz5Kw4YNadKkCc8991zu/vft28fw4cOJiIigd+/ebNmy5Zjj5x1hTktL44477iAuLo7IyEgGDBhAWlracaUXSUlJDB8+nKioKNq0acPTTz+du7+pU6cycuRIrrzySsLDw4mPj2fNmjUluhbFbduyZUuWLFkCQFZWFg8++CCtW7cmPDycHj16kJiYeNw5LVy4kG7duhEREUFsbCxTp04tcd+IiIhI6TjneH75Ns77+2fs+PUwT47uzhV9Wvg7rAqnhLu07rkHDh8+dtnhw97yMrBixQrS09OLHCl+4IEHWLlyJd988w1r167liy++YPr06bntv/zyCwcPHmTnzp08++yz3HLLLfz6668A3HLLLYSEhPDzzz8ze/ZsZs+eXehxJk2axJdffsny5cvZv38/Dz/8MAEBx//KjBo1iubNm5OUlMR//vMf7r77bj766KPc9nfffZfLL7+cAwcOMHz4cG699dYSX4+SbvvYY48xb948Fi1axKFDh5g9ezZ16x7/VKrQ0FBeeOEFDhw4wMKFC3nyySd5++23SxyPiIiIlNxjH25kyrvr6ds6mvcnnMbQTk38HZJfKOEure3bS7e8lPbt20dMTAxBQYVPIDN37lwmT55Mw4YNadCgAVOmTOHFF1/Mba9VqxaTJ0+mVq1anHvuuYSFhbFhwwaysrJ44403uP/++wkNDaVjx46MGzeuwGNkZ2cze/ZsnnjiCZo1a0ZgYCD9+vUjODj4mPUSExNZtmwZM2bMICQkhK5du3LttdceE8+AAQM499xzCQwMZOzYsaxdu7bE16Ok2z7zzDNMnz6ddu3aYWZ06dKF6Ojo49YbOHAgnTp1IiAggM6dOzNq1CiWLl1a4nhERESkeOlHvYkVrujTggdGdOS5q3rRMCLEz1H5jxLu0mpRyNcghS0vpejoaPbu3VvkTBlJSUnExcXlvo+LiyMpKemYfeRN2OvWrUtKSgp79uwhMzOT2NjYY7YtyN69e0lPT6d169ZFxpuUlERUVBTh4eHH7HPnzp257xs3bnxMLOnp6SWeCaSk2yYmJhYbK8CqVasYNGgQDRo0IDIyklmzZrF3794SxSIiIiJFO3j4KLe/8jXXvbCG7GxHk8g6jO4Th1n1nF+7pJRwl9YDD0D+UoW6db3lZaBv376EhIQUWebQtGlTfvrpp9z327dvp2nTpsXuu0GDBgQFBeXWNudsW5CYmBhCQkKOq/EuKJb9+/eTnJx8zD6bNWtWbDxlKTY2tthYAa644gqGDx9OYmIiBw8e5MYbb8Q5VwERioiIVG9LN+5h8N+WsnDdz/SMi0L/u/6PEu7SGj0annoK4uLAzPvzqae85WUgMjKS+++/n1tuuYW3336bw4cPc/ToUd577z3uvPNOwKuZnj59Onv27GHv3r3cf//9jBkzpth9BwYGctFFFzF16lQOHz7Md999x/PPP1/gugEBAYwfP56JEyeSlJREVlYWK1asyL05M0dsbCz9+vXjrrvuIj09nXXr1vHss88yuoyuR0lde+213HfffWzatAnnHOvWrWPfvn3HrZecnExUVBQhISF88cUXvPzyyxUap4iISHVz+Egm97z1X8bN/oKIkFq8dXN/bj/rZAIDavaodl560uSJGD26zBLsgkycOJFGjRoxffp0Ro8enTvrxj2+GzPvvfdeDh06ROfOnQG49NJLuffee0u073/+859cffXVNG7cmPbt23P11VfzySefFLjuI488wl133UWvXr1ISUmhS5cuvP/++8etN2/ePG688UaaNm1K/fr1+fOf/8zZZ599gmd/YiZOnEhGRgaDBw9m7969tG/fnrfeeuu49WbOnMkdd9zBrbfeyumnn87IkSM5cOBAhcYqIiJSnWRmOxI27OG637XijsHtCKkVWPxGNYxV96/Te/bs6Qqahu7777+nQ4cO5Xrs5OTkY2qbpfpRH9cM+fu5Iv79kIqXkJBQoodpSdWlPi476UezeGHFNsb1a0lwUCCpGZmEBleOcVx/9rOZfemc65l/eeW4MiIiIiJSJaxNPMAdr69l8+4U4qJDOSe+caVJtisrXR0RERERKdaRzGz+/tEmnly6hQZhwTw/vjent23g77CqBCXcIiIiIlKsP76xjre+3sklPZpz3/mnEFmnlr9DqjIqfJYSM5ttZrvN7NtC2i8ws3Vm9o2ZrTGzAXnaxpnZJt+r4Ce2iIiIiEiZOJqVTWqG9/yL6087iWfH9eSRS7so2S4lf0wLOAcYUkT7R0AX51xXYDzwDICZRQFTgD5Ab2CKmdX/LYFU9xtGRaTsZWdn+zsEEZEK8V3SIS745+dMfmc9AB2aRHBmh0Z+jqpqqvCE2zn3KbC/iPYU979MOBRy500/B/jQObffOfcr8CFFJ+5FCgkJYd++fUq6RaREnHMcOXKEnTt3Ehoa6u9wRETKzdGsbP62ZCPD/7mM3ckZDI5Xkv1bVcoabjMbATwENATO8y1uBiTmWW2Hb9kJad68OTt27GDPnj0nHGdx0tPTCQkJKbf9i/+pj2uGnH4OCgoiMjKSmJgYf4ckIlIuNu9O4ffzvua7nw8xolszpgw7hXp1a/s7rCqvUibczrm3gLfM7DRgGnAWUNDjigocnjaz64HrARo1akRCQkI5RVq0lJQUwsLC/HJsqRjq45ohbz/v3LnTz9FIeUlJSfHb/xdSMdTHxdufns3+Q+nc3j2Ybg0P8M0Xy/0dUqlVxn6ulAl3Dufcp2bW2sxi8Ea0B+Zpbg4kFLLdU8BT4D34xl+Tn2uC/epPfVwzqJ9rBvVz9ac+LtjaxAO8/mUi0y7oiJlx4WBHQBV+LHtl7Gd/3DRZJDNrY2bm+7k7UBvYB7wPDDaz+r6bJQf7lomIiIhIKaUfzeKh975nxMzPWfLdbpIOpgNU6WS7sqrwEW4zm4c3Uh1jZjvwZh6pBeCcmwVcDFxpZkeBNOAy302U+81sGrDat6v7nXOF3nwpIiIiIgVbvW0/f/zPOrbuTWVU71juOrcDESGa6q+8VHjC7ZwbVUz7DGBGIW2zgdnlEZeIiIhITXA0K5s/vPoNAC9d04cBJ+tG8PJWqWu4RURERKRsLN+ylx5x9QkOCuSZcT2JrV+X0GClghWh0tVwi4iIiEjZ+TX1CBNf+4Yrnl7Fiyt+AqB94wgl2xVIV1pERESkGnLOsei/vzDl3W85cPgotw5qw5hT4/wdVo2khFtERESkGvrr+xuYmbCFTs0ieWF8H05pGuHvkGosJdwiIiIi1UR2tiMjM5s6tQM5v3NTIuvU4poBrQgKVBWxPynhFhEREakGNu9O5q43/0uLqFAeHdmFU5pGaFS7klDCLSIiIlKFHcnM5smELfzrk83UqR3IyJ6x/g5J8lHCLSIiIlJFbfglmdvmfcXGXSkM69KUyeefQoPwYH+HJfko4RYRERGpoiLr1MIwZl/VkzPaN/J3OFIIVdCLiIiIVBHOOd7778/8ft7XOOdoHBnC4gm/U7JdyWmEW0RERKQK2HkgjSnvfMuS73cT3zSC/alHiA4Lxsz8HZoUQwm3iIiISCWWmZXNnOXbeOzDjTgH95zbgav7t9RUf1WIEm4RERGRSuxIVjbPfb6N3q2imHZBR2Kj6vo7JCklJdwiIiIilUxy+lGeXfYjN57emrq1g3j7lv7EhNVW+UgVVeKE28zqAb8DegONgRBgP7AR+Nw59025RCgiIiJSQzjnWPztL0ydv57dyRl0aV6PQe0baqq/Kq7YhNvMTgNuA84HagPbgb1ABtANuBoINbMfgWeBfznnDpVbxCIiIiLVUOL+w0x+51s+2bCHU5pE8O+xPekaW8/fYUkZKDLhNrMlQHfgDWAEsDx/Mm3edxvtgKHApcAkMxvnnFtQPiGLiIiIVD//95+1/HfHQe49rwNX9dNNkdVJcSPci4ELnXMpha3gnHPAD77X42bWG4gquxBFREREqqeVW/dxcsMwosOCeWBEJ+rWDqRJZB1/hyVlrMiE2zn3SGl36Jz74sTDEREREan+9qVk8OCiH3jjqx1c97tW3HPeKbRuEObvsKSclPi7CjP72MzaF9LW1sw+LruwRERERKqfrGzHSyt/YtAjCby7die3DGrNxLPb+TssKWelmRZwIBBRSFsEcNpvjkZERESkGnvsww3865Mt9D0pmmkXxtOmYbi/Q5IKUNp5uF3+BWZWGzgD+KVMIhIRERGpRg4ePkrKkUya1avD2FNb0rZROMO7NNWc2jVIkSUlZjbFzLLMLAsv2V6Z8z7P8jTgIeClCohXREREpEpwzvHGlzs487EE/u/1tQA0jgzhgq7NlGzXMMWNcC/Cm3PbgL8DjwLb8q1zBPjBOfdZmUcnIiIiUgV9//MhJr/zLau3/UrX2HrcfW4Hf4ckflTcLCWrgdUAZpYMLHDO7auIwERERESqoiXf7eKGl74kIiSIGRd34tIesQQEaES7JivuwTfmm2cb59zzJd1p3u1EREREqjvnHHuSM2gYEULf1tFc97uTuPH0k6hXt7a/Q5NKoLhpATea2bVmFlqSnZlZDzN7AfjTbw9NREREpPL7/udDXPbvlVz21EoyMrMIDQ7iT0PbK9mWXMXVcP8JmAo8YWYfAMuBb/HqujOAekAroAcwBIgFngFml1O8IiIiIpXCwbSjPP7hRl5c+RMRIUH83zntqRWgx7HL8Yqr4X7DzN4EzgKuBCYATfjf9ICGd9Pkl8C/gRedc7vLL1wRERER/9u6J4WR/17BvtQjjO7TgkmD22lEWwpV7DzcvlrsD30vzKwJ0BgIAfYDPzrnjpRnkCIiIiKVwaH0o0SE1CIuOpTB8Y25oncLOjaL9HdYUsmV9sE3OOd+Bn4uh1hEREREKqX9qUf46/sbeH/9LyyZeDpRobV5cEQnf4clVUSJC43M7H4zKzBBN7NoM3u97MISERER8b/MrGyeX76NgX/9hNfWJDKiWzNqBWqKPymd0oxw3w6cb2ZXOue+zVloZiOAWUByWQcnIiIi4i8pGZlc8uRyfvglmf5topk6LJ6TG4X7OyypgkpzK20X4ACwxsz+ZGYxZjYXeAN409cuIiIiUqWlZmQCEBYcRN/W0cwa052XrumjZFtOWIkTbufcNufcGcCdwBRgJzAAOMc5d5NzLrWcYhQREREpd+lHs/j7R5s49aGP2Lw7BYApw+IZ0rEJZiojkRNXqpsmzSwM6AwEA7vwZioJK+U+ZgPnA7udcx0LaB8N/NH3NgW4yTm31te2Da90JQvIdM71LM2xRURERPJzzrH421+YvvB7dh5I49xOjalbO9DfYUk1UuKE28wGAc/5tjkf+BiYAfzHzOYBtznnfi3BruYA/wReKKT9R+B059yvZjYUeArok6d9kHNub0njFhERESlMdrbjqjmr+XTjHto3Dufl6/rQr3WMv8OSaqY0I9xLgFeAW/Mk1rf7HowzG/gv0Ly4nTjnPjWzlkW0L8/zdmVJ9ikiIiJSGikZmYQFBxEQYPSKq8/ZHRoyqncLggL1pEgpe6X5rbrMOTc6/yi2c24p3g2TC8o0Ms81wHt5Dwd8YGZfmtn15XA8ERERqcaOZmUz5/Mf6f+Xj1m2yfvC/LYzT2Zs35ZKtqXcmPcgyQo+qDfCvaCgGu486wwCZgIDnHP7fMuaOueSzKwh3pMvb3POfVrAttcD1wM0atSoxyuvvFL2J1ECKSkphIWVqsRdqhj1cc2gfq4Z1M/V3xeJKby9LYCkVEd8dACjOwTTNExJdnXjz7/LgwYN+rKgewxL/aTJimBmnYFngKE5yTaAcy7J9+duM3sL6A0cl3A7557Cq/2mZ8+ebuDAgRUR9nESEhLw17GlYqiPawb1c82gfq7e7nhtLW+sTyUuug5PX3wKZ3VoqJlHqqnK+He50iXcZtYCb17vsc65jXmWhwIBzrlk38+Dgfv9FKaIiIhUcgfTjhIWHERggNGnVRSBqbuYNvY0goM0A4lUrApPuH0zmgwEYsxsB96c3rUAnHOzgMlANDDT98kzZ/q/RsBbvmVBwMvOucUVHb+IiIhUbplZ2cxbncjjH27kznPacXnvFozsFUvD1C1KtsUvKjzhds6NKqb9WuDaApZvRU+zFBERkSJ8tmkP0xZ8x8ZdKfRpFUWn5pH+DknkxBJu84aZm+A9vCazbEMSERERKb0/z1/Pc59vo0VUXWaN6cE58Y1Upy2VQmmfNHkuXglIVyAQ76bFr8zsKeBT59xLZR+iiIiISMF+TT1C7aAAQoODGNiuIY0iQri6f0uVjkilUuK5cMzsSuBd4Ae8KffybrsJb85sERERkXJ3JDObZz7byul//YSZCZsBOL1tA248vbWSbal0SjPCfQ/wV+fcXWYWiPeY9xzrgUllGpmIiIhIPs45PvxuFw+99wM/7k3ltLYNuKBrM3+HJVKk0iTccXgPmylIOhDx28MRERERKdyMxRuYtXQLbRqG8dzVvRjUrqG/QxIpVmkS7kSgG/BxAW09gc1lEpGIiIhIHj8fTCPAjEYRIVzYrSnN6tdhVK9YPYpdqozS/KY+C0wxszFAHd8yM7MzgTuBp8s6OBEREam5UjMyeeyDDQx6JIEZ7/0AQPvGEYw9NU7JtlQppRnhngHEAs8DWb5ly/FmK/m3c+7vZRybiIiI1ECZWdm8/uUOHvtwI3uSMxjWpSl/OLutv8MSOWElTridcw64xcweA87CexrkfuDjvI9gFxEREfkt/v7xZv7+0SZ6xNXn32N70L1FfX+HJPKblPrBN865LcCWcohFREREaqj1SQcJMKNDE69kpEPjcIZ0bKwH10i1UGTCbWanlWZnzrlPf1s4IiIiUpP8fDCNR97fyJtf72BQu4bMvqoXDcKDGdqpib9DEykzxY1wJwAOyPl46fK0Wb734NVzi4iIiBQpOf0oTyZs4dllP+IcXPe7k7hlUBt/hyVSLopLuDvl+bkJMBtYDLwJ7AYaAhcD5wDjyyNAERERqX5eXZ3IzIQtXNC1KZMGtyM2qq6/QxIpN0Um3M659Tk/m9mDwAvOuXvzrbbYzKYDE4AlZR+iiIiIVHXOORZ/+wvBtQI4o30jxpwaR59W0XRqHunv0ETKXWkmsTwTWFpI21Jg4G+ORkRERKqd1dv2c/GTy7lp7le8tHI7ACG1ApVsS41RmllK9gMXUPDj3Uf42kVEREQA2LInhb+89wMffreLhuHB/OWiTlzSo7m/wxKpcKVJuP8C/NPMWgLv8r8a7guAocCtZR2ciIiIVF3f/3yIFVv2MWlwW8YPaEXd2qWejVikWijNg29mmtlO4G7gn75tM4FvgIucc2+XT4giIiJSFSSnH+XpT7cSUacW1/7uJM7r1IR+rWOICq3t79BE/KpUHzWdc+8A75hZANAA2OOcyy6XyERERKRKyMjM4uVV2/nHx5vZn3qEUb1jATAzJdsinMCTJgF8SfauMo5FREREqphlm/Zy11vrSNyfRr/W0fxpaHs6N6/n77BEKhUVU4mIiEipOOc4kpVNcFAgocGBhAfX4oXxnfjdyTF6FLtIAZRwi4iISImtTTzAjMU/EBddl4cu6ky3FvVZ+PsBSrRFiqCEW0RERIq1dU8Kj3ywgUX//YWo0NqcE984t03JtkjRlHCLiIhIkd74cgd3vrGO4KAAfn/myVz3u1aEh9Tyd1giVYYSbhERETnOwcNHOZR+lNiouvQ5KYqxp8Zxy6A2NAgP9ndoIlVOqRJuM7sMuA5oC4Tkb3fONSyjuERERMQP0o5kMWf5Np5M2EzXFvV5YXxvmtevy9Th8f4OTaTKKnHCbWZXALOBOcAZvp8DgOHAAeCFcohPREREKsDRrGxeXZ3I3z/axO7kDM5o35BJg9v5OyyRaqE0I9z/B0zDe8T79cBM59xXZhYOfAgcLof4REREpAI8v3wb0xd+T8+4+vzziu70bhXl75BEqo3SJNwnA58757LMLAuIAHDOJZvZDOBx4JFyiFFERETKmHOOTzbsJiQokH5tYrisVywnNQhlULuGmnVEpIyVJuE+COTcKbET6AAk+N4bEF12YYmIiEh5WbV1H399fwNrfvqVM9s3pF+bGMJDanFG+0b+Dk2kWipNwr0G6Ay8D7wLTDazTOAIMBlYVfbhiYiISFn5LukQMxb/wNKNe2gUEcwDIzoysmesv8MSqfZKk3A/BMT5fp7s+3kmEAisxqvrFhERkUpqfdJB1u44wF1D2zOuX0tCagX6OySRGqHECbdzbiWw0vfzAeACMwsGgp1zh8opPhERETlBifsP88RHm4hvGsHV/VsxolszzunYmAg9tEakQpVmWsBwIMw593POMudcBpBhZk2AZOdcSjnEKCIiIqWw+1A6//xkM/O+2I6ZERdVF4CgwAAiAgP8HJ1IzVOakpJn8W6cvK6AtqlAJHB5GcQkIiIiJ2juqp+4f/53ZGU7RvaK5bYz2tAkso6/wxKp0UrzMfc0YGEhbYt87cUys9lmttvMvi2kfbSZrfO9lptZlzxtQ8xsg5ltNrM/lSJ2ERGRautQ+lEOHj4KQKuYUM7r1ISP7jidB0d0UrItUgmUJuGOpPCH26QD9Uu4nznAkCLafwROd851xnvQzlMAZhYI/AsYCpwCjDKzU0p4TBERkWrn8JFMZiZs5nczPuGJjzYB0K91DI9d1pW46FA/RyciOUpTUrIJOA/4oIC2c4EtJdmJc+5TM2tZRPvyPG9XAs19P/cGNjvntgKY2SvABcB3JTmuiIhIdZF+NIu5q7bzZMJm9qYcYVC7BlzUvZm/wxKRQpQm4f4HMMvMjuCNUv8MNAHGAbcAN5V5dHAN8J7v52ZAYp62HUCfgjYys+vxTVPYqFEjEhISyiG04qWkpPjt2FIx1Mc1g/q5ZqhK/fz8+gw+ScykQ1QAN/QJ4eT6h9m76WsSNvk7ssqtKvWxnLjK2M+lmRbwaTNrBNwFTMzTlA7c65x7uiwDM7NBeAn3gJxFBYVVSKxP4StF6dmzpxs4cGBZhlZiCQkJ+OvYUjHUxzWD+rlmqMz9nJmVzZtf76RbbD1ObhROy46pXHcwjX6tY/wdWpVSmftYyk5l7OfSjHDjnJtuZv8A+uI9yn0fsMI5d7AsgzKzzsAzwFDn3D7f4h1A3sdhNQeSyvK4IiIilUlWtmP+2iSe+GgTP+5N5YbTTuKuczvQMiaUljGq0RapKkqVcAP4kuvF5RALAGbWAngTGOuc25inaTVwspm1AnbiTUF4RXnFISIi4k/vr/+Fv76/gc27U2jfOJx/j+3B4FMa+TssETkBRSbcZnYusMw5d8j3c5Gcc4uKW8fM5gEDgRgz2wFMAWr5tp+F99j4aGCmmQFkOud6OucyzexW4H28x8nPds6tL+54IiIiVYVzXqWkmfHVT79iwMzR3RkS35iAgIIqK0WkKihuhHsBcCrwhe/noji8RLjolZwbVUz7tcC1hbQtwpvzW0REpNpwzvHR97t5fMlG/u+cdgxs15AJZ7XlziHtCVSiLVLlFZdwt8KbjSTnZxERESkjzjkSNu7h8Q83sm7HQVpE1cU3yE2d2sWOYYlIFVFkwu2c+wnAzIKBMcAC59zaighMRESkurvppa9YvP4Xmtevw8MXd2ZE92bUCizNM+lEpCoo0U2TzrkMM7sHWFbO8YiIiFRbzjlWbN1Hz7goagcFMKRjY05v14CLuzendpASbZHqqjSzlKwCegBLyykWERGRask5x4ot+/jbkk18sW0/f72kM5f2jOXCbno6pEhNUJqE+07gZd+TJhcBu8j34Bnn3OEyjE1ERKTKW7FlH48v2cgXP+6nUUQw918Qz/CuTf0dlohUoNKOcAP8HXiikHV0h4eIiIiPc45pC75jX2oGfx4ez2W9Ygmppf8qRWqa0iTc4ynkUeoiIiLiWbFlH09/tpXHR3Ylsm4tnhzTnUYRIUq0RWqwEifczrk55RiHiIhIlebVaG9k1Y/7aRgezJa9KXRvUZ+4aD2CXaSmK/Wj3c3sFLybJ2Pxnvb4i5m1AXY555LLOkAREZHKLP1oFuNmf5GbaE8ZdgqjerfQiLaI5Co04TazYOdcRp73YcBs4BLgqG/bxcAvwIPAdmBSuUYrIiJSCTjn2LInhTYNwwmpFUirmFCGdGysRFtEClTUpJ8rzCwuz/vHgH7AmUA4kPdZs4uAIWUfnoiISOXhnOPTjXu4ZNYKzvnbZ/y0LxWAv1zcmav7t1KyLSIFKqqk5GPgCzMb7ZxbAlwE3O6c+8TM8v+L8hMQd9weREREqoGcR7D//aNNfL39AE0jQ5g6PJ7GkSH+Dk1EqoBCE27n3CQzext4BmgP1AH2FbJ6OJBV9uGJiIj43+7kDK5/YQ0Nw0N4YERHLunRnOAgjWaLSMkUedOkc26ZmXX3vV0NXIlXt53fJcDyMo5NRETEL7KzHR9+v4sVW/YxdXg8jSJCePm6U+nSvJ4ewS4ipVbsLCV5nh55L7DEzJYAr+PNyX2umf0BL+E+rdyiFBERqQDZ2Y73vv2Ff3y8iR9+SaZldF0OHm5LZN1a9GoZ5e/wRKSKKs083MvM7EzgL8A/8W6a/DOwEjjLObe6fEIUEREpfzuTsznnb5+yaXcKJzUI5fHLujCsc9P/b+/O46uqzv2Pf57MTAFCmBMIo8ggyAxKCTjhVIcWh2Idbr3a69BRa2uttrW97VV776/WWmtbSrEqztapDkXiiIDgiBMzCVOAyBBCQobn98c5aBozHTjn7CTn+369zivn7LVz9hdWSB5W1l6LlGSNaIvI4YloHW53fw2YZmbtgK7Arloj4CIi+M1z9QAAIABJREFUIq1KZXUNW3eXk5vVnqx2Rtf2afzu/KM5ZVRvkpOs6TcQEWmGiDe+AXD3/cD+KGcRERGJi4qqah5dsYk7C1aTnpLM89/5Eu1SjAe/OSXoaCLSBkVUcJvZROAsoC9Qdy0kd/dzoxVMREQk2sorq1mwdCN/fHktW3aXMzq3C9+aORjTYLaIxFCzC+7wzZG/AbYBa4EDsQolIiISC8+t3MpPn/yAiQOyuOWrR3Hs4GxM1baIxFgkI9zfB34LfM/dPUZ5REREomb3/krmv76erh3SuGByf04d1Zu+XdoxXiuOiEgcRVJwpwNPq9gWEZGWbmdpBXNfW8f81zewt6KKc8bnAP1JSU5SsS0icRdJwT2P0Pbu/4pNFBERkcN335KN3PzUB5RXVXPKyN5cMWMQI/p0DjqWiCSwSAru64A7whvfvAjsqtPu7v6HqCUTERFppsKSMjJSk+neKZ287PacPKoXV+QPZnCPjkFHExGJqOCeCcwBOoWf1+WACm4REYmb1cV7uXPRGv7xzmYumpLHjacPZ+qgbKYOyg46mojIZyIpuO8ElgDfBla7e2VsIomIiDTu/U27+f2i1Ty7cisZKclcPDWP/5w2MOhYIiL1iqTg7gNc4e4fxiqMiIhIc/z5lbW8unoHV80YzCXHDCCrQ1rQkUREGpQUwbn/AkbHKoiIiEh93J1FHxdzzl2LWbl5NwDXn3Ikr/1wJt8/8QgV2yLS4kUywn07cJeZtaP+myZx9w+iFUxERBJbdY3z7Ptb+f2i1XywZQ99OmewfW8FAD0y6252LCLSckVScB9cDvDnwM/qtBmhmyaToxFKREQSm7tz1p2v8W7RbgZmd+CWrx7FmWP6kpYSyS9mRURahkgK7hkxSyEiIgmv7EAVT7+7ha+Oy8HMmD0uh8u/NIhZI3uRnKTt10Wk9Wp2we3uL8UyiIiIJKZdZQf42+sbmPf6Oj4tq2Rg9w6M65/F16fkBR1NRCQqml1wm1n7ps5x97LDiyMiIolib3klty9cxX1LNrLvQDXHDevBFTMGMa6/tl4XkbYlkiklpYTmaTdGc7hFRKRR+w9U0y4tmbSUJJ55byvHD+/Jf+UPYlivzKCjiYjERCQF93/wxYI7CzgRGA7cHK1QIiLS9ry/aTd/KFjD24W7WHRNPukpySz8/nQyUjVWIyJtWyRzuOc10PR/ZnYnMKI572Nmc4HTgGJ3H1lP+zDgr8BY4MfufluttvXAXqAaqHL38c3NLyIi8efuLF6zkz+8tIZXVu2gU3oKX5/Sn8rqGtJSklRsi0hCiGSEuzGPAg8A323GufOAO4D5DbSXAN8CzmygfYa774g0oIiIxN+SdSV87c9LyO6Yzg9mHcEFk/uTmZEadCwRkbiKVsE9Aahozonu/rKZ5TXSXgwUm9mp0YkmIiLxUlFVzeNvbaLsQDWXHDOASQOy+O15YzhpRC+NZotIwjL3pu6DDJ9odks9h9OAI4HjgP/n7tc0873ygKfqm1JS65yfAqV1ppSsAz4lNJf8j+5+dwOfexlwGUDPnj3HLViwoDmxoq60tJSOHTsGcm2JD/VxYlA/N21/lbOosJLn11exq8IZlpXEdRMyMGs962ern9s+9XFiCLKfZ8yYsby+Kc+RjHDPrudYOVBEaApIvcVvlB3j7pvNrAfwgpl95O4v1z0pXIjfDTB+/HjPz8+PQ7QvKigoIKhrS3yojxOD+rlxT76zmesfe4+95VUcM7gb35w+iGMHZ7eqYhvUz4lAfZwYWmI/R3LT5IBYBmlmhs3hj8Vm9hgwEfhCwS0iIrG1bsc+UpKM3Kz2DMjuwLQh2Xxz+iCOyukSdDQRkRYnKegAzWVmHcys08HnhJYjfD/YVCIiieXtwl18857lzPxNAbcvXAXAyL6duXPOOBXbIiINaHSE28wujOTN3L2hlUdqv+f9QD6QbWZFwE1Aavjz7zKzXsCbQCZQY2bfIbTOdzbwWPhXlCnAfe7+bCT5RETk0Lyyajt3vLiaJetKyMxI4Yr8QVw0NS/oWCIirUJTU0rmNeM9at912WTB7e7nN9G+Fcipp2kPMLoZeUREJAoqq2tISTLMjIUfFrOxpIwbTj2S8yb2o2N6tBa5EhFp+5r6jtmpifaJwA3ADOCjqCQSEZFAlVZUsWDpRua+uo7bzhnN1EHZfO/EoVx/ypGkpbSamYgiIi1GowW3u++r77iZTQN+DJwAvAOcCzwc9XQiIhI3xXvLmffaeu55YwN7y6uYNCDrs7WztVmNiMihi+h3gmZ2IqFCexqwBDjD3Z+KRTAREYmfmhrnrN+/zpbd+5k1sheXfWkQY3J1E6SISDQ0q+A2sy8TKrQnAC8BJ7j7wlgGExGR2HF3lq3/lEeWF/HLs0aSkpzEL88aSV63DuRldwg6nohIm9LUKiXnAtcDI4FngWPd/fV4BBMRkeirrnGeW7mVP768lncKd9G1fSrfmDaAoT07kX9Ej6DjiYi0SU2NcN9PaBWSAmAlcKaZndnAue7u10Uxm4iIRNGmXfv52p/eYMPOMvpltefmM0bw1XG5tEtLDjqaiEib1lTBvZFQwT0w/GiMAyq4RURakB2lFXy8dS/HDM6md2YGo3O6cN2sYZw0ohfJSa1r63URkdaqqVVK8uKUQ0REomjt9lL+9Mo6HllRRIe0ZJZcfzxpKUncfv7RQUcTEUk42rlARKQN+WjrHm577hMWfrSN1OQkvjI2h0unDdD62SIiAVLBLSLSylVV11BWWU1mRiplB6pZvqGEq2cM5utT8ujeKT3oeCIiCU8Ft4hIK1VaUcWDywqZ+9o6pg3J5ldnH8XYfl1Z/KPjPtuwRkREgqeCW0Skldmyez/zXlvPfUs3sre8igl5XTn+yJ6ftavYFhFpWVRwi4i0MncuWsO9SzZw8sjeXDptAEf36xp0JBERaYQKbhGRFqymxnnpk+386ZW1fPeEoUzIy+LKGYO57EsDyc1qH3Q8ERFphogKbjPrA5wG5AAZdZq18Y2ISJSUV1bz6IpN/OXVtazZvo9emRmU7DsAQK/Odb/9iohIS9bsgtvMziK082QyUAwcqHOKNr4REYkCd+e0373K6uJSRvbN5LfnjeGUUb1JTdbSfiIirVEkI9z/DTwPXOzuJTHKIyKSkD7eupdHVxTxg1nDSE4yvnXcEHp0SmfSgCzMtCOkiEhrFknBnQtcrWJbRCQ63EPzs//y6jpeWbWDjNQkzhjTl+F9Mvny6D5BxxMRkSiJpOB+HTgC+FeMsoiIJIxNu/Zz8dylrCoupUendK496Qi+NrEfXTukBR1NRESiLJKC+3vAvWZWCrwA7Kp7gruXRSuYiEhbU7y3nNXbSpk6OJtemRkMyO7Af+UP4rSj+mjrdRGRNiySgvvd8Me/ErpBsj7abUFEpI4PNu/hL6+u48l3NpPZLpXFP5pJanISd184PuhoIiISB5EU3P9Bw4W2iIjU8XbhLv7nnx+xeO1O2qUmc97EXC45ZoBWGxERSTDNLrjdfV4Mc4iItAn7KqqoqKohq0MaFZXVrN+5jx+ePIzzJ/Sjc/vUoOOJiEgAIt5pMrz5zRQgCygBFrv75mgHExFpTTbt2s/819dz/9KNnHl0X35+xkgmDsji5R/M0Ii2iEiCa7DgNjNzd6/1Ohn4HfCf/Ptc7Wozu5vQkoE1MUsqItICrdj4KX95dR3Pvr8VgJNH9uKso/sCYGakJmsNbRGRRNfYCPczZnaeu+8Ov/4ZoXnc1wMPANuAnsC5wM+BncCNMcwqItIiVFXXkBIetX5wWSGvfLKdS6cN4MIpefTt0i7gdCIi0tI0VnDvA940szPdfSVwIXCDu99W65yNwK1m5sC3UMEtIm1Yyb4D3L90I/MXr+euC8ZxdL+uXHPSEdx4+nDap0U8Q09ERBJEgz8h3P2rZvZ14J9AP6AHny8NWNe74XYRkTanaG8NP3zkXR57axMVVTVMG5JNSlJohDu7Y3rA6UREpKVrdEjG3e8xsxfDLz8BzgOer+fU84CPo5xNRCRwFVXV/Grpfqp8E2eP7cvFUwdwRK9OQccSEZFWpMnfgbr7pvDTXwALzKwf8DChOdw9gNnADEJFt4hIq7a3vJKHlxfx0ifbmXvRBNJTkrlqTAbnzZqmbddFROSQRLIO94NmtovQzZO/BVKBSmA5MMvdX4hNRBGR2Fu/Yx9/W7yeh94sorSiirH9urBz3wG6d0rnyG7JKrZFROSQRXSXj7s/DzxvZklANrBDSwGKSGu3ZO1OzvvTG6QkGaeO6s3FxwxgTG6XoGOJiEgbcUi31YeL7OIoZxERiYuyA1U8/tZmkpPg3An9GNe/K9eedARfHZtDj8yMoOOJiEgb02jBbWa3ALe7e1H4eWPc3a+LXjQRkegq+rSMexZvYMGyQnbvr2TmsB6cO6EfKclJXJE/OOh4IiLSRjU1wj0buBcoCj9vjAMquEWkRfpDwRpufe4jzIxZI3px8TF5jO/fNehYIiKSAJpaFnBAfc8Ph5nNBU4Dit19ZD3tw4C/AmOBH9feaMfMZhG6YTMZ+LO7/zoamUSk7dl/oJp/vL2JYwZnk5vVnrH9unD59EF8fXJ/+mg3SBERiaMgtkabB9wBzG+gvYTQrpVn1j5oZsnA74ETCI24LzOzJ9z9g9hFFZHWpujTMu55YwMPLCtkV1klP5h1BFfkD2bSwG5MGtgt6HgiIpKAml1wm9kvgWx3v7yetruA7e7+k6bex91fNrO8RtqLgWIzO7VO00RgtbuvDV9zAXAGoIJbRHB3vr3gbZ56dzMAJ43oxcVT85g4ICvgZCIikugiGeE+H7ixgbZXgJ8DTRbch6EvUFjrdREwqb4Tzewy4DKAnj17UlBQEMNYDSstLQ3s2hIf6uNgVVQ77++oZlzP0Ley8l0VnJyXysx+KXRrt5f9G9/jpY2Hfx31c2JQP7d96uPE0BL7OZKCuw+wqYG2zeH2WLJ6jnl9J7r73cDdAOPHj/f8/PwYxmpYQUEBQV1b4kN9HIzCks+njezeX8nz353E0J6diFVXqJ8Tg/q57VMfJ4aW2M+RFNxbCd3IuKietrHA9qgkalgRkFvrdQ6hQl9EEsSmXfu56R8rWfjRNpLMmDUyNG1kSI+OQUcTERFpUCQF94PAjWb2kbs/ffCgmZ1CaCrJ3dEOV8cyYIiZDSA00n4e8LUYX1NEAlZaUcXmXfsZ2rMTndulsqp4L1fmD2bO5H707qzVRkREpOWLpOC+ERgDPGlmO4EtQG8gC3ieZs7fNrP7gXwg28yKgJuAVAB3v8vMegFvAplAjZl9Bxju7nvM7CrgOULLAs5195UR5BeRVmTt9lLmL97Aw8uL6NulHc9+Zxod01MouCYfs/pmmImIiLRMzS643b0cONHMTgJmAN2AncBCd38hgvc5v4n2rYSmi9TX9gzwTHOvJSKtz9J1JfzuxVW8smoHqcnGaUf14aKpeZ8V2Sq2RUSktYl4HW53f47QKLOISFR8uu8AKclGp4xUNu/az6ptpXzvhKGcP7Ef3TulBx1PRETksESyDvfwps7RJjQiEon3N+1m/uL1/OPtzXz3hKF8c/ogTjuqN6ce1ZvU5KSg44mIiERFJCPc79PAMny1JB9GFhFJEP94exN/e309Kzbuol1qMl8Zl8PMYT0ASFGhLSIibUwkBfeMeo5lASeGH9+OSiIRaZN2l1XSuX0qAI+s2MSnZZX85LThfHVcDp3bpQacTkREJHYiuWnypQaaHjOzXwDnAE9FJZWItAnuzuK1O7ln8QYWflTMomvy6dulHb89dwyd26WSlKQbIEVEpO2L+KbJBiwCHo3Se4lIK7evoopHVhRxz+INrCoupUv7VC6ZmkdaeLpI1w5pAScUERGJn2gV3KcCu6L0XiLSSpVXVpORmsze8ip+9uQHjOiTya1fPYrTR/chI1W3eIiISGKKZJWSB+s5nAYMA4YA10crlIi0HpXVNTy/chvzF68nJdm499LJ9OqcwYvfn07/bh2CjiciIhK4SEa4u9dzrBx4BfheeFMaEUkQ2/aUc9+Sjdy/dCPFeyvI6dqOCyb3x90xMxXbIiIiYZHcNFnfKiUikkDcnRqH5CTjyXc2c/uLq5g+tDu/mtyf/CN6kKybIEVERL4gWnO4RaQN21NeyWMrNnHPGxu4/EsDmT0+l3Mm5HLC8J4ayRYREWlCowW3md0SwXu5u193mHlEpAX5YPMe/r5kA4+/tYmyA9WMzu1CdsfQVuuZGalkZmj9bBERkaY0NcI9O4L3ckAFt0grV1Pjn62Pfe3D77C6uJQzxvThgsn9OSqnS8DpREREWp9GC253HxCvICISrMKSMu5dspGn39vM09+aRmZGKr85ZzS9MjPo0l7rZouIiBwqzeEWSWDVNU7Bx8Xcu2Qjiz4uxoAThvdkz/5KMjNSGdYrM+iIIiIirV5EBbeZDQSuBY4FsoASQssC3ubua6MfT0Ri4eDSfWu3l/KNv71J907pXD1zCOdPzKV353ZBxxMREWlTItn4ZhyhLdzLgaeAbUBP4CvAHDOb4e4rYpJSRA6bu7N0XQl/X7KR9JQkbps9miE9O/H3b0xi0sAsUsPbrouIiEh0RTLCfRvwFnCyu5cdPGhm7YFnwu0zoxtPRA7XnvJKHl1exL1LNrKquJTMjBTOn9Tvs/Zjh2QHmE5ERKTti6TgngicU7vYBnD3MjO7DXggqslE5LAcnDZyV8Ea7ixYw+icztzylaM4fXQf2qUlBx1PREQkYURScO8HujXQlkVoqomIBGj/gWqefHcz976xge8cP5QZw3pw4ZQ8Th7Zm1E5nYOOJyIikpAiKbifBn5tZmvd/dWDB83sWOBXwJPRDicizbNq217uXbKRR1YUsbe8isE9OlLjDkCvzhn06pwRcEIREZHEFUnB/T3gH8BLZrad0E2TPcKP14HvRz+eiDTk4JSRmhrn4r8uY/veCk4e1Ys5k/ozIa8rZhZ0RBERESGCgtvddwLHmtksYALQG9gCLHH352OUT0Tq2LBzH/ct3chLH2/nyauPJTU5idvPP5q8bu3pFt52XURERFqOiDe+cfdngWdjkEVEGlBZXcPCD4u5d8kGXlm1g+Qk4/gje7CrrJLundIZ179r0BFFRESkAYe002R4KcBvAMOArcB8d98QzWAi8vm0kTfXf8o3/76c3p0z+O7xQzl3Qq7mZYuIiLQSjRbcZvYb4HR3H1rrWCdgGTAE+BToDHzfzCa6+yexDCuSCA5ut37fko0MyO7ADacNZ/LALP72HxM5ZlA3UrRBjYiISKvS1Aj3DODvdY5dAwwFLnX3uWbWHXgB+Anw9ehHFEkM2/aU88CyQhYs3cjm3eV075TOhAFZAJgZ04d2DzihiIiIHIqmCu48YHmdY18BPnD3uQDuvj08Ev6z6McTadtqapykpNBqIr95/mMefLOIaUOyufH04Rx3ZE9tty4iItIGNFVwp1BrQxszywKOBH5f57z1QK+oJhNpw4oPjmYvK+SuC8YxKqczV88cwpUzBtO/W4eg44mIiEgUNVVwfwLkAwvDr08Lf3yuznk9gJLoxRJpe2pqnFdW7+C+JRv414fFVNc4xwzu9tkGNblZ7QNOKCIiIrHQVMF9B/AnM+tMaKObbwHrgLrrbp8IvB/9eCKtX2V1DanJSVRU1XDVfStITU7i0mMHcN7EfgzI1mi2iIhIW9dowe3u88ysN3Al0AVYAVzp7pUHzwnfNHkGmsMt8pnao9kbdpbxz29Po11aMvddOpmhvTqSnpIcdEQRERGJkybX4Xb3XwG/aqR9O5q/LQJA8d5yHgzPzS76dD9ZHdKYPS6HiqoaMlKTGZXTOeiIIiIiEmeHtPGNiHyuusaprA4V1EvWlnDb858wdVA3fjBrGCeN6KnRbBERkQSnglvkEG3dXc4/Vh/gx28s4muT+nHljMGcOKIni67J19xsERER+UzcC24zm0totZNidx9ZT7sBvwVOAcqAi919RbitGngvfOpGd/9yfFKLfG7RR8Xcu2QDL35UTI3DMYMzGdEnE4D0lGQV2yIiIvJvghjhnkdo9ZP5DbSfTGjb+CHAJOAP4Y8A+919TKwDitS1s7SCbh3TAbjnjQ28W7Sby6cPYkDNZs45ZXLA6URERKQli3vB7e4vm1leI6ecAcx3dwfeMLMuZtbb3bfEJaBIWGV1DQs/3Mb9Swt5bfUOFl2TT25We379lVF0bZ9GanISBQVbg44pIiIiLZx5eNONuF40VHA/1cCUkqeAX7v7q+HXC4Hr3P1NM6sC3gaqwuc83sD7XwZcBtCzZ89xCxYsiMmfoymlpaV07NgxkGvLodtzwHl2XSWvbqpizwEnK8OY1jeF4/qnkplm/3au+jgxqJ8Tg/q57VMfJ4Yg+3nGjBnL3X183eMt8aZJq+fYwf8V9HP3zWY2EHjRzN5z9zVfONn9buBugPHjx3t+fn7MwjamoKCAoK4tkamoqmb73gpyurZn+94Krn35RaYN6cHXJuUyfWgPkpPq+7JUHycK9XNiUD+3ferjxNAS+7klFtxFQG6t1znAZgB3P/hxrZkVAEcDXyi4RZpr1ba9LFhWyKMrihjasxMPXD6F7p3SWXbD8WRmpAYdT0RERNqAllhwPwFcZWYLCN0sudvdt5hZV6DM3SvMLBs4BrglyKDSei38cBt3Fqxh+YZPSU02Thzei/Mmfv7/PBXbIiIiEi1BLAt4P5APZJtZEXATkArg7ncBzxBaEnA1oWUBLwl/6pHAH82sBkgiNIf7g/iml9bK3Xl/0x4Gdu9Ah/QUCkvK+LTsANefMoyzx+aQHV6BRERERCTaglil5Pwm2h24sp7jrwOjYpVL2qbd+yt54u1N3L+0kA+27OHXZ4/ivIn9mDO5PxdNzSO07LuIiIhI7LTEKSUih62iqpofPfIeT7+3hYqqGkb0yeTmM0dy8qjeAKQmJwWcUERERBKFCm5pM4r3lvNe0W6OO7In6SnJFO+tYPb4HM4d349ROZ2DjiciIiIJSgW3tGpV1TW8vGo7C5YWsvCjYlKTjeU3nECH9BTu+cZETRkRERGRwKngllbrpU+2c93D77J1TznZHdP4xrEDOHdCLh3SQ1/WKrZFRESkJVDBLa1GeWU1z63cSl63DozO7ULfLu0Y1rsTP/3yCI47sofmZYuIiEiLpIJbWryVm3fz4LJCHn97M7v3V/L1yf0ZnduFwT06Mu+SiUHHExEREWmUCm5p0S792zL+9WExaSlJzBrRi3Mn5DJlYLegY4mIiIg0mwpuaTFqapwl60p45r0t3HT6cFKSk8g/ogfThnTnjDF96NI+LeiIIiIiIhFTwS2B27J7P48sL+LBN4vYWFJGp4wULpjcnyN6deKCyf2DjiciIiJyWFRwS6De37SbL9/xKjUOUwZ247snDGHWiN60S0sOOpqIiIhIVKjglrj6eOteHnyzkK7tU7lq5hCG987k+ycewWlH9aZ/tw5BxxMRERGJOhXcEnN7yit58p3NPPhmEe8U7iI12ThnfC4ASUnGlTMGB5xQREREJHZUcEtMuPtnG8/c/OQHPLS8iKE9O3LDqUdy1tF96dYxPeCEIiIiIvGhgluiavOu0A2QDy0v4q4LxjG8TyaXTx/InMn9GZ3TWbs/ioiISMJRwS2H7UBVDc+t3MpDy4t4ZdV2PHwD5IHqGgAG9+gUcEIRERGR4KjglkPi7uzeX0mX9mlUVtfww0fepXO7VK6eMZjZ43PJzWofdEQRERGRFkEFt0SkZN8BHn9rEw8tL8Ld+ee3p9EhPYXHrzyGgd07kpykKSMiIiIitanglmZZsfFT/vTyWv714TYqq51RfTtzzvhcahySDYb01LQRERERkfqo4JYGrdleSlb7NLp2SGP9jn0sXVfChVPymD0+h2G9MoOOJyIiItIqqOCWf7O3vJKn3t3CQ28WsmLjLn508jAunz6I047qw2lH9SEtJSnoiCIiIiKtigpuAaCmxrn24Xd5+r3NlFfWMLhHR3508jDOGtsXQIW2iIiIyCFSwZ3ACkvKWLa+hLPH5pCUZJRXVnP22Bxmj8thTG4XrZktIiIiEgUquBNM2YEq/vneVh5aXsgba0tISTJmDutBl/Zp/H7O2KDjiYiIiLQ5KrgTyIsfbePq+95i34Fq+ndrz/dPGMrZ43Lo0j4t6GgiIiIibZYK7jZs0679PLaiiBF9OjNjWA+G9+7MKaN6M3t8LhPyumrKiIiIiEgcqOBuY/YfqOa5lVt5eHkRr63ZgTtcPn0gM4b1oFfnDG6dPTroiCIiIiIJRQV3G/O1P7/BWxt3kZvVjm8fN4SvjM3RNusiIiIiAVLB3YodnDLy7MqtPHT5VNqlJfPt44aQnpLMpAFZJGmbdREREZHAqeBuZcoOVPHs+1t5ZEURr6/ZiTtMGdiNHaUV5Ga1J/+IHkFHFBEREZFaVHC3AjU1TlllNR3TU1hTvI/vPfiOpoyIiIiItBIquFuwwpIyHllRxCMripg8oBu3zh7NyL6ZPHrFVMbkdNGUEREREZFWQAV3C/T0u1v42+L1LF1XghlMHdTts6kiZsbYfl2DDSgiIiIizaaCuwWornGWriv57EbHNzeUsH1vBdeedARnHt2Xvl3aBR1RRERERA6RCu4ArdleyiPLi3jsrU1s2V3O/f85mSmDunHdrGHceNpwbUwjIiIi0gao4A7Apl37ufLeFbxduIskgy8N7c71pxzJ0f26AJCRmhxwQhERERGJlqR4X9DM5ppZsZm930C7mdntZrbazN41s7G12i4ys1Xhx0XxSx2he++FvDymz5wJeXlU3fN3/vXBNp58ZzMAPTql0y41metPGcYbPzqOeZdM5PTRfVRoi4iIiLRBQYxwzwPuAOY30H4yMCT8mAT8AZhkZlnATcB4wIHlZvaEu38a88SRuPdeuOz3XIBDAAAKu0lEQVQyKCvDADZsoPIbl/LESVex4aQzOX10H1KTk7j/sslBJxURERGROIj7CLe7vwyUNHLKGcB8D3kD6GJmvYGTgBfcvSRcZL8AzIp94gj9+MdQVvZvh9pVVvA/yx/g4f+aGlAoEREREQlKS5zD3RcorPW6KHysoeNfYGaXAZcB9OzZk4KCgpgErc/0jRup71bHjK2beOmVl+OWQ+KjtLQ0rl9fEgz1c2JQP7d96uPE0BL7uSUW3PXVq97I8S8edL8buBtg/Pjxnp+fH7VwTerXDzZs+MJh69ePuOaQuCgoKFC/JgD1c2JQP7d96uPE0BL7Oe5TSpqhCMit9ToH2NzI8Zbll7+E9nW2Wm/fPnRcRERERBJOSyy4nwAuDK9WMhnY7e5bgOeAE82sq5l1BU4MH2tZ5syBu++G/v1xM+jfP/R6zpygk4mIiIhIAOI+pcTM7gfygWwzKyK08kgqgLvfBTwDnAKsBsqAS8JtJWZ2M7As/FY/d/fGbr4Mzpw5MGcOL7XAX2mIiIiISHzFveB29/ObaHfgygba5gJzY5FLRERERCQWWuKUEhERERGRNkMFt4iIiIhIDKngFhERERGJIRXcIiIiIiIxpIJbRERERCSGVHCLiIiIiMSQCm4RERERkRhSwS0iIiIiEkMquEVEREREYshCGzu2XWa2HdgQ0OWzgR0BXVviQ32cGNTPiUH93PapjxNDkP3c39271z3Y5gvuIJnZm+4+PugcEjvq48Sgfk4M6ue2T32cGFpiP2tKiYiIiIhIDKngFhERERGJIRXcsXV30AEk5tTHiUH9nBjUz22f+jgxtLh+1hxuEREREZEY0gi3iIiIiEgMqeA+BGY218yKzez9BtrNzG43s9Vm9q6Zja3VdpGZrQo/LopfaonEofaxmY0xs8VmtjJ8/Nz4JpdIHM6/5XB7ppltMrM74pNYDsVhfs/uZ2bPm9mHZvaBmeXFK7dE5jD7+Zbw9+0Pw+dY/JJLczWjj4eFfwZXmNk1ddpmmdnH4f7/YXwSf04F96GZB8xqpP1kYEj4cRnwBwAzywJuAiYBE4GbzKxrTJPKoZrHIfQxUAZc6O4jwp///8ysSwxzyuGZx6H180E3Ay/FJJlE0zwOvZ/nA7e6+5GEvm8XxyijHL55HNrP5qnAMcBRwEhgAjA9lkHlkM2j8T4uAb4F3Fb7oJklA78n9DUwHDjfzIbHKGO9VHAfAnd/mVCnNuQMYL6HvAF0MbPewEnAC+5e4u6fAi/Q+BeOBORQ+9jdP3H3VeH32Ezoh/MXFsCXluEw/i1jZuOAnsDzsU8qh+NQ+zn8AznF3V8Iv0+pu5fFIbIcgsP49+xABpAGpAOpwLZY55XINdXH7l7s7suAyjpNE4HV7r7W3Q8ACwh9PcSNCu7Y6AsU1npdFD7W0HFpfZrsSzObSOgb+Jo45pLoqrefzSwJ+A1wbSCpJNoa+vc8FNhlZo+a2Vtmdmt4pExap3r72d0XA4uALeHHc+7+YQD5JHYCr79UcMdGfXO/vJHj0vo02pfhUZN7gEvcvSZuqSTaGurnK4Bn3L2wnnZpfRrq5xRgGnANoWkGA4GL4xdLoqzefjazwcCRQA6hImymmX0prskk1gKvv1Rwx0YRkFvrdQ6wuZHj0vo02Jdmlgk8DdwQ/rWltF4N9fMU4CozW09oruCFZvbr+MeTKGnse/Zb4V9DVwGPA2Pr+XxpHRrq57OAN8JThkqBfwKTA8gnsRN4/aWCOzaeIPQD2MxsMrDb3bcAzwEnmlnX8M2SJ4aPSetTbx+bWRrwGKF5gg8FG1GioN5+dvc57t7P3fMIjX7Od/e43/UuUdPQ9+xlQFczO3gfxkzgg6BCymFrqJ83AtPNLMXMUgndMKkpJW3LMmCImQ0I/5w+j9DXQ9ykxPNibYWZ3Q/kA9lmVkRo5ZFUAHe/C3gGOAVYTWjVikvCbSVmdjOhjgf4ubs3doOHBORQ+xg4B/gS0M3MLg4fu9jd345beGm2w+hnaUUO43t2dXhpsYXhZeKWA3+K+x9AmuUw/j0/TOg/U+8RmmbwrLs/Gdfw0ixN9bGZ9QLeBDKBGjP7DjDc3feY2VWEBjmTgbnuvjKu2bXTpIiIiIhI7GhKiYiIiIhIDKngFhERERGJIRXcIiIiIiIxpIJbRERERCSGVHCLiIiIiMSQCm4RkYCYmTfjkW9mF4efdww6c7SE10J+x8wuqnXsB2aWX+e8dmZWbGbT4h5SRCRKtCygiEhAwptvHNQOeBH4BaGdSg/6AEgHBgFL3b0mfgljx8zOBW4FBrl7ZfjYDuAOd/9pnXN/Ahzn7vnxzikiEg3a+EZEJCDu/sbB57VGr9fUPl7L9vikiptvAfccLLabMA/4mZmNcvf3YhtLRCT6NKVERKSFq29KiZllmNktZlZoZhXh6Rmn1Pm89WZ2m5n90My2mNluM/tNeDrHKWa20sz2mtnjZta11uflh693opk9ZWb7zGyjmX2znmznmNl74QyFZvZLM2t0MMfMBgNTCe3w91lWoBtwU+3pNADuXkhoh94LD+XvT0QkaCq4RURap4eBi4H/Bk4nVJA+YWZj6px3HjCR0DbWtwDfA/4XuBn4CfBNYDrwq3qu8RfgXeBs4J/AH8zstIONZnYi8ACwAjgD+B1wDXBHE9mPA/YB79Q6dhawO3zNKeHHilrtrwPHN/G+IiItkqaUiIi0MmZ2HHAqkO/uL4UPP29mQ4EfA7NrnV4OzHb3auBZMzsDuBoY4u7rwu83GriIUPFd2z/d/frw8+fMbCBwA/BU+NjPgQJ3P3jj47NmBvArM/uFuxc18EcYB3xYez66u79lZlVAUQNTat4BrjazDHcvb+jvRkSkJdIIt4hI63M8sBV4zcxSDj6AhcD4OucWhIvtg1YD6w8W27WOdTeztDqf+1id148C48ws2cySgbHAQ3XOeYDQz5YpjeTvBexopL0+O4BkoHuEnyciEjiNcIuItD7ZhIrW+m44rK7zeled1wcaOGZAWvj5QcV1zism9HMjO/w6FdhW55yDr7PqCx6WAZQ10l6filqfKyLSqqjgFhFpfUqATcCZMb5Oj3peV/H56HRlPef0DH8saeR9Swj9hyESXZrxviIiLZKmlIiItD4LCRWspe7+Zt1HFK9zVj2vl7t7dXiaynL+fb44wDlADbC4kff9GBhQz/EDNDyCnQfsdPedTYUWEWlpNMItItL6vAA8B7xgZv8DrAQygTFAhrv/KErXOdnMfgm8RGilkhMIrUZy0E2Ebqb8K7AAGEVo9ZM/NXLDJMBrwI1m1t3da68v/hFwqpk9C5QCH7v73nDbeEIrlYiItDoa4RYRaWU8tEXw2cBc4DuEiu8/ErpR8dUoXupSQjdGPg6cBlzp7k/UyvE8oWUHxwNPhrP8BriqifctIDQ1ZFad49cSWi7waULLHI4DCN8QehzwyGH9aUREAqKt3UVE5N+EN5xZBIxy9/djdI3fAoPd/dRmnHsS8CDQx933xSKPiEgsaYRbRESCcCuQH147vCnfBf5PxbaItFYquEVEJO7Cc7y/AfRu7Dwza0foBsz/jUcuEZFY0JQSEREREZEY0gi3iIiIiEgMqeAWEREREYkhFdwiIiIiIjGkgltEREREJIZUcIuIiIiIxJAKbhERERGRGPr/Tmp67yTEmBsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "f = lambda t,x: np.tan(x)+1\n",
    "x0=1\n",
    "t0=1\n",
    "tf=1.1\n",
    "n=100\n",
    "t,x=RK38(f,x0,t0,tf,n)\n",
    "plt.figure(figsize=(12,6))\n",
    "plt.plot(t,x,'--',label=\"Solución Numérica con RK38, n = \"+str(n))\n",
    "plt.plot(t[0],x[0],'ro', label=\"Condición Inicial\")\n",
    "plt.xlabel(\"Tiempo (t)\",fontsize=15)\n",
    "plt.ylabel(\"Solución Numérica de x(t)\",fontsize=15)\n",
    "plt.title(\"Método de Runge-Kutta (RK4)\",fontsize=15)\n",
    "plt.legend(fontsize=12)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora tenemos **muchos métodos numéricos a nuestra disposición** para resolver muchas ecuaciones diferenciales.\n",
    "\n",
    "¡Buenas noticias!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tarea"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construya y programe los siguientes métodos numéricos usando la Tabla del Carnicero correspondiente:\n",
    "\n",
    "- Heun.\n",
    "- Punto Medio.\n",
    "- Ralston.\n",
    "- Dormand-Prince (Opcional).\n",
    "\n",
    "Una todas esas definiciones y los métodos ya construídos (Euler, RK4, RK38) en un módulo llamado ```Solver_EDO```.\n",
    "\n",
    "Resuelva las siguientes ecuaciones, usando cada método numérico  en el módulo creado y comparando los errores correspondientes.\n",
    "\n",
    "1.\n",
    "$$\\frac{dP}{dt}=kP,\\text{   }P(0)=2,\\text{   en  }T=[0,100]$$\n",
    "2.\n",
    "$$\\frac{dP}{dt}=kP\\biggl(1-\\frac{P}{S}\\biggr),\\text{   }P(0)=2,\\text{   en  }T=[0,100],\\text{   }S=100000$$\n",
    "3.\n",
    "$$\\frac{dv}{dt}=g-\\frac{c_d}{m}v^{2},\\text{   }v(0)=-3,\\text{   en  }T=[0,30]$$\n",
    "4.\n",
    "$$\\frac{dx}{dt}=e^{-t^2}\\sinh(x), \\text{  }x(-5)=0$$\n",
    "5.\n",
    "$$\\frac{dx}{dt}=\\tanh(xt)+\\ln(t^2)e^{-x^2}, \\text{  }x(-1)=2$$\n",
    "\n",
    "Decida cuál es el método numérico más apropiado para cada ecuación."
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}