Diferencias finitas.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Diferencias finitas\n",
    "\n",
    "\n",
    "Una forma de resolver el problema de Laplace numéricamente es en términos del llamado método de diferencias finitas. En este, el dominio es discretizado introduciendo una grilla rectangular de puntos equiespaciados, separados por una distancia $\\varepsilon$ suficientemente pequeña como para que el potencial que esperamos obtener varíe suavemente entre puntos vecinos de la grilla.\n",
    "\n",
    "Aproximamos entonces el valor de la derivada direccional en un punto,  a lo largo de uno de los ejes de la grilla, como \n",
    "$$\n",
    "\\frac{\\partial \\phi(\\vec{x}_0)}{\\partial \\zeta_i} \\approx \\frac{\\phi(\\vec{x}_0+\\varepsilon\\hat{e}_i) -\\phi(\\vec{x}_0) }{\\varepsilon}\n",
    "$$\n",
    "Notemos que en esta expresión, la derivada direccional en un sentido y en el opuesto no tienen por qué coincidir. Una mejor aproximación consistiría en asignar este valor a la derivada direccional *en el punto medio* entre $\\vec{x}_i$ y $\\vec{x}_i + \\varepsilon \\hat{e}_i$.\n",
    "\n",
    "Bajo esta hipótesis, podemos aplicar nuevamente la misma regla para obtener una estimación de la derivada segunda en el punto $x_0$:\n",
    "$$\n",
    "\\frac{\\partial^2 \\phi(\\vec{x}_0)}{\\partial \\zeta_i^2} \\approx \\frac{\n",
    "\\frac{\\partial \\phi(\\vec{x}_0+\\varepsilon/2 \\check{e}_i)}{\\partial \\zeta_i} -\n",
    "\\frac{\\partial \\phi(\\vec{x}_0-\\varepsilon/2 \\check{e}_i)}{\\partial \\zeta_i}\n",
    "}{\\varepsilon}\\approx \\frac{\\phi(\\vec{x}_0+\\varepsilon \\check{e}_i)+\\phi(\\vec{x}_0-\\varepsilon \\check{e}_i)\n",
    "-2\\phi(\\vec{x}_0)}{\\varepsilon^2}\n",
    "$$\n",
    "\n",
    "El laplaciano de $\\phi(\\vec{x})$ en $x_0$ se aproxima entonces como\n",
    "$$\n",
    "\\nabla^2 \\phi(\\vec{x}_0)= \\frac{\\overline{\\phi(\\vec{x}_0)}-\\phi(\\vec{x}_0)}{\\varepsilon^2/n} \n",
    "$$\n",
    "donde $n$ es el número de vecinos y $\\overline{\\phi(\\vec{x}_0)}=\\frac{\\sum_{i,\\pm}\\phi(\\vec{x}_0\\pm\\varepsilon\\check{e}_{i})}{n}$ es el *valor promedio* de $\\varphi$ sobre todos los puntos vecinos a $x_0$\n",
    "\n",
    "De esta manera, una solución para esta discretización de la ecuación de Laplace debe satisfacer que en cada punto, el valor de $\\varphi$ coincide con el promedio sobre sus vecinos. \n",
    "Una condición semejante podría haberse obtenido partiendo de la propiedad de la Ec. de Laplace continua, que establece que el promedio de $\\phi(\\vec{x})$ sobre una esfera coincide con el valor que toma en su centro.\n",
    "\n",
    "\n",
    "Las condiciones de borde tipo Dirichlet se imponen ahora fijando el valor del potencial sobre los puntos de la grilla cercanos a la superficie que la define. De la misma manera, la Ec. de Poisson toma la forma discretizada\n",
    "$$\n",
    "\\phi(x_0)= \\overline{\\phi(x_0)} + \\frac{\\varepsilon^2}{n \\epsilon_0} \\rho(x_0)\n",
    "$$\n",
    "\n",
    "con lo que podemos construir la función de Green del problema remplazando $\\rho(x)$ por una carga $Q=\\varepsilon^3 \\rho(x)$\n",
    "\n",
    "\n",
    "Las ecuaciones discretizadas definen un problema lineal con restricciones, que se puede resolver en principio por inversión directa. Sin embargo, en la práctica se suele preferir el uso del llamado \"método de relajación autoconsistente\": se parte de una aproximación a la solución $\\phi_0$, que cumple con las condiciones de borde.\n",
    "A partir de esta aproximación, se produce una nueva aproximación vía el remplazo\n",
    "$$\n",
    "\\phi_1(\\vec{x})=\\overline{\\phi_0}(\\vec{x}) + \\frac{\\varepsilon^2}{n \\epsilon_0} \\rho(x_0)\n",
    "$$\n",
    "\n",
    "Iterando este procedimiento hasta que $\\max_{\\vec{x}}|\\phi_{n-1}(\\vec{x})-\\phi_{n}(\\vec{x})|$ sea menor que una cierta tolerancia, encontramos la solución buscada.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Algunas parametrizaciones de curvas donde definir las condiciones de frontera"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [],
   "source": [
    "def curva_frontera_semicirc1(t):\n",
    "    t = t * np.pi\n",
    "    y = int(n/2. * (1.-.5*np.cos(t)))\n",
    "    x = int(n/2. * (1.+.5*np.sin(t)))\n",
    "    return (x,y)\n",
    "\n",
    "def curva_frontera_semicirc2(t):\n",
    "    t = t * np.pi\n",
    "    y = int(n/2. * (1.-.5*np.cos(t)))\n",
    "    x = int(n/2. * (1.-.5*np.sin(t)))\n",
    "    return (x,y)\n",
    "\n",
    "\n",
    "\n",
    "def curva_frontera_quart_circ1(t):\n",
    "    t = t * np.pi\n",
    "    return (int(0.25 * n * np.cos(t)), int(0.25 * n * np.sin(t)))\n",
    "\n",
    "def curva_frontera_quart_circ2(t):\n",
    "    t = t * np.pi\n",
    "    return (int(n-0.25 * n * np.cos(t)), int(n-0.25 * n * np.sin(t)))\n",
    "\n",
    "\n",
    "def curva_frontera_rect1(t):\n",
    "    t = t * np.pi\n",
    "    if t< .25* 3.1416:        \n",
    "        x = .25\n",
    "        y = .5 + .25*np.tan(t)\n",
    "    elif t> .75* 3.1416:\n",
    "        x = .75\n",
    "        y = .5 - .25*np.tan(t)\n",
    "    else:\n",
    "        x = .5 - .25/np.tan(t)\n",
    "        y = .75\n",
    "    return (int(n*x),int(n*y))\n",
    "\n",
    "\n",
    "def curva_frontera_rect2(t):\n",
    "    t = t * np.pi\n",
    "    if t< .25* 3.1416:        \n",
    "        x = .25\n",
    "        y = .5 - .25*np.tan(t)\n",
    "    elif t> .75* 3.1416:\n",
    "        x = .75\n",
    "        y = .5 + .25*np.tan(t)\n",
    "    else:\n",
    "        x = .5 - .25/np.tan(t)\n",
    "        y = .25\n",
    "    return (int(n*x),int(n*y))\n",
    "\n",
    "\n",
    "def curva_frontera_tri1(t):\n",
    "    x = 0.25+.5*t\n",
    "    y = .5 * (1- 1./np.sqrt(12.)) \n",
    "    return (int(n*x),int(n*y))\n",
    "\n",
    "def curva_frontera_tri2(t):\n",
    "    x = 0.25+.5*t\n",
    "    t = np.abs(.5-t)\n",
    "    y = .5 * (1- 1./np.sqrt(12.))  + np.sqrt(.75)*(.5-t)\n",
    "    return (int(n*x),int(n*y))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parámetros"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "lado=10 # tamaño físico del sistema (en cm, por ejemplo)\n",
    "n=50    # número de puntos  en cada dirección \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuI5Wle3/H39/ldzq1O3asv29PsRddNZEE0Q2JcEHEVjIqbPxJQUDZGGAKJrq5B1kD0r4B/iOgfITB4BZcNYV1QJBhlowkhsmR2FZzd2dkdZ3dm+lLdVV3XU+f2u3zzx+90V3X3OdXVVTVdXac+L2i6q/r8Tj1TU/Xppz7n+T2PuTsiIjJdwlkPQERETp/CXURkCincRUSmkMJdRGQKKdxFRKaQwl1EZAo9MdzN7HfN7K6ZvXrgfYtm9hdm9vXR7wvv7jBFRORpHGXm/vvADz3yvk8Bn3f3DwKfH70tIiLPCTvKTUxm9j7gT939w6O3Xwe+z91vm9lV4K/c/UPv5kBFROTo4mNed9ndbwOMAv7SpAea2UvASwCWpP8oXZn4UBERGWNw68a6u688zTXHDfcjc/eXgZcB6teu+3v/zSff7Q8pIjJVvvYrn3zraa857mqZO6M6htHvd4/5PCIi8i44brj/CfDx0Z8/Dvzx6QxHREROw1GWQn4G+GvgQ2Z2w8x+Bvg14AfN7OvAD47eFhGR58QTO3d3/4kJf/XRUx6LiIicEt2hKiIyhRTuIiJTSOEuIjKFFO4iIlNI4S4iMoUU7iIiU0jhLiIyhRTuIiJTSOEuIjKFFO4iIlNI4S4iMoUU7iIiU0jhLiIyhRTuIiJTSOEuIjKFFO4iIlNI4S4iMoUU7iIiU0jhLiIyhRTuIiJTSOEuIjKFFO4iIlNI4S4iMoUU7iIiU0jhLiIyhRTuIiJTSOEuIjKFFO4iIlNI4S4iMoUU7iIiU0jhLiIyhRTuIiJTSOEuIjKFFO4iIlNI4S4iMoVOFO5m9gtm9mUze9XMPmNm9dMamIiIHN+xw93MrgE/B7zo7h8GIuDHT2tgIiJyfCetZWKgYWYx0ARunXxIIiJyUscOd3e/Cfw68DZwG9h29z9/9HFm9pKZvWJmrxR7e8cfqYiIHNlJapkF4GPA+4H3AC0z+8lHH+fuL7v7i+7+YtRqHX+kIiJyZCepZX4A+Ia7r7l7BnwO+J7TGZaIiJzEScL9beC7zaxpZgZ8FHjtdIYlIiIncZLO/QvAZ4EvAX83eq6XT2lcIiJyAvFJLnb3XwV+9ZTGIueIZWDF5L/3BDx6duMRkYedKNzlAioh6kGUHf4w7ztlzSjqgD2TkYnIAQp3ORqHkEHUc6yE2pZT2/HxDzXoLxgZRsicvGm4vtJEnil9y8mTlRB3IeQQhtBcK6mtD0g3B2Mf7sGI9+oMFhJ6K4G445SpUTTQLF7kGVG4yxPFexAKsMKpb5XEuznJzhBf24LNrccvCIH67iWS5RniXoP+UmDYrmb+edPw5Nn/N4hcNAp3eaKiAT5wIofupYjuciC+nlLfaNH6ch37+5uENIUwWnzlJeUbbxFtLtDsLBN3WwwWE7orgWQPiqR6Tu1JKvLuUbjLE3kMRWwUTcAdj0r6cyWdodFducz81TbpqzcI/Ryr1zAzrF6n7HTg9V1q3feQ7MwQ9Zv0FwLDOSPkkDcMT8/6v05kOinc5VCOQ4AyKfGFnLQ5ZGmxw/X2Bt085fWVq9xZbjC78q20v7JB9I3bhHodiyKi2VnK4ZDyG28TNudo7V0iudRksJzSXQmYQ5lBrlm8yKlTuMtEZezk7ZJ4Zkg6O+Bbl9Z5b2uD72q/xVzcIy8D/7f9Lfyv2W9jY26O4cIy84stkr97C7Z3qG5cBgsRfm8TdjqknSvEu7OEXpPBQsRgwUgyKBpGmaIXXEVOicJdJioThxKS2SGtdMBCo8/1+j3m4h6Xwi5rzPDBxl3emL/E1y7X2Rs2CMMGC9l14tVtfDB88FwPMjsYcTen+XaHKGuRdKNqFg+Uw0NufDIoamiGL3JECnc5nIEFJ4lLAOohB2Ap6tInYSXdZbHRozYzpLuQsHc9Iek3aSYBxi+DJwxLSGOSbkk0dChj8pbRXzDIJw8lDNXTixyVwl2OrWEZsZWkIWe52WF1KWLQj9i9bnjUGH+RQ9xz0t2CMCgJudNYzyg6gXQnPKhyHlWkqKcXeQoKdzm2ulXT7OVkl3w2sNqZJbQz+mWKFUYoJhTopdNch2gQUbuXEQ0Kks0ulBOm+kDRjKuefjEwmA9VT980ygT19CJjKNzl2GoHwv32cIFmMqS/6HQ7CYMliHqPp645RN1qvXzcLXFPqK8VRKXjWzv41vbYjxVdWmaml5PsNUm6Kd3lUU8fQ1GvtjwYK6DwlwtJ4S7HFkalemTV73EoedCqBKoj0x/hQN42wgAcI3/BKBoNsnZCoxYRFmbHfzAHjyOigVO/l0Mx6ukXjZBPTm83V08vF5LCXZ49g7IOZWLEXeitGMN2QtloEyZsI+w4yV5BGPh+T78XiLsRRX1yuA/mDXMoMt0VKxeLwl3OTgT5TLUKxoHdFyYnr5XQuGdEA6jdywglxHsFyWafMh1/nRskO7UHPX1QTy8XiMJdzpZBWatm8dGAicsnLfeHevp0J6e2uoeVTtjeGf9irBnJQptkt0XSre2vp48hb6JZvEw1hbs8H8KoNpnEH+7pG/cSirRNY7WLbRt+5y4UJVY7cJurO35nnbRzhaQzS9RrMViI6C8ayQ54PHn6XtTQ7pVyrinc5XwY19PPJJT1GWozCclsC7+xinf7RHOzWFx9aXuWU9xYxTa3ae5dIVlpknRr9JYCxYQXWT1AyE27V8q5pnCX8+XRnv56RNZqUGsnNJoptrZNcXOVUKtjrSaWxESLC3i3S/n6myS7l4k780TdBp6M3+vArXohdjAfCPlo3xv19HLOKNzl/Lnf08dG3IP+UjWLLxoz1Nop6UwTv3WXcmOD0GiOrglYrYa/cwu2tmnuXsbS8b2LB0iWG/s9vaunl/NH4S7nVwR5q5rFR0DnWkTWqlWz+FZKWNuhXN8AH73YGoCZJrhT3lqFbh/K8vHnjWPSF64Q784Rek2GB3p67V4p54XCXc63A6tt4h4MFiKymUA5msUnC20mLcGxXobfXIVOlzDbftDTU5QUN1YJm9u0uldIRz19dyV6sHtl3mTsTVoizwuFu0yHUM3ibQhxDzpXA8Nmjdp8On6S7WDDgkar6unLgz19FO339F99k2TnMvHufDWLX4zpLRrJLhR1o6yhWbw8lxTuMlU8hSw2oh4M543h/OTHhkGgqB/o6W/eodzYIJqdxZIEa7WwWo1ida3q6feukuw2ibt1esv3Z/FO0bLJ+9CLnBGFu0yfAEULyhzskP3hI6DzwuM9fXFjFbKMBxvluEOW4V97k2TvCnFnjtBvkc1Xh4XbrlPWjKKOZvHy3FC4y9TyuPo1SZmOevrFUU9fn6HWTqi1W5BN+FchBGxYUl8fEGU1oj70F40MIwx9/x+ER8cSRnvQa4Yvz4jCXS6uUL0wahnEBp33BIatOlk7JeqP38HMHEK/JGDEg5IoK/GQkHQKhm2rHjBG3jCSXD29PDsKd7nY7PGePmsbUX/8gnYrId12ap2CqFtWG5rdzSgjSHbDxNAu07Df02dO0VRPL+8uhbsI7Pf0GUR9o5y074w72YxRbBpxx6lvZcRbA9JBjufF+HXzQDlTI/RaZAvq6eXZULiLHOAJ5IdtGFZWe9v0F40w4xAl1CgJawW+twed7uPXmBHCEq23S4b9JqGfMrjf02dOUZuc7mWMeno5FoW7yNMYracPGUQHe/pWQq2dEPLxM3fLSjxNCDk0NjIICcleQW85YOMvAapDStTTy3Eo3EWelkGZVnvbPOjpZ4x8pjnxJClySLuP9PQ1I+rFE2fmbtBfCurp5VgU7iLHdaCnj7vQvXxI8rpTf9DT54SsJN7OSO/s4fGEk6SCEe01GC7E9JZNPb08lROFu5nNA78NfJhqA49/7e5/fRoDEzkvPIFs1qobpiacJBUN7UBPH5Nu59S3MsgL/Bs3xl4XokBzd4V0uUXcq9NfCmQYdn8Wr6mZHOKkXx6/BfyZu/8LM0uB5imMSeT8scNPbsqTh3v6tJWQ19o07nSJVpYo37kFWU5oNfdvhMoKyq99g3hriVZnibjXYrAQ01sOhM7Ef0fAqpOk1NNfbMcOdzObBb4X+FcA7j4EhqczLJEpM66nb1m1B/1sSq1Vh1trlBtbhNk2oVarLmvUKTd3YHuXevcqyfIMcbfBYC5M7N+LumFu6ukvuJPM3D8ArAG/Z2bfAXwR+IS77x18kJm9BLwEEM8tnODDiUyBR3r6vauBrFUjbyXUWwnR+izlO7fx/oDQnql2qJyfo+z3Kf/+baKteZp7yySLjYlnwBa1iOF8TG+p6um12uZiOkm4x8B3AT/r7l8ws98CPgX8x4MPcveXgZcB6teuT/xJUuQiud/TRz0YzlWrbYpGi3Q2pdasweo9irWNh06LsijC796D7Q61y0sQRRAefzHW04R06f7ulQHDqt0r1dNfKCf5X30DuOHuXxi9/VmqcBeRozAomlVdE3dh70rEsGUUrZj6TEo0Pwfbu49eAl7i/SG+tQ3DDIKBHQx5J15ZptVZfKint45W21wkxw53d181s3fM7EPu/jrwUeArpzc0kYvBY0b72YBjZK2IvN4knUtJtievUYh3l+DWGr6xRWjPPOjpvSgoN7dhe4d69z0kyy3iXoPeaLVNyJy8YYe+ACzn30l/SPtZ4NOjlTJvAj998iGJXEAGRaM6LjDqOt3LEcO2kS5M+BZ1iLuNI/T0bz3o6aNui+FCQnclkJRQpNXH1Cx+Op0o3N39b4EXT2ksIheex5C3jTCoZvF5c/wNTgBW+MM9/e17FPc2CM1mdZJUFBHaM5Qb27DTob53lWRlhqjfpL8YGM4aIYO8qVn8NNLLKyLPG4OyXh0mMnHfGa/Oin3Q0zdHPf3CLOXNVdjZefjxeY5/8wZxfpWZbk7IWiR7cTWLdygSDu/i7ZC/k+eSwl3keRWqE5wmyeIDPX0zIm80qc2m1Js1QjZ5N7LSjHS7IBpAyGKGs8Zg3oiyyel9/yQpzfDPD4W7yHk1rqefMfLmLKEYv+rYgaRTEDInypz6RkbcCyR7EXljwhGBVDddJaWppz9HFO4i59zjPb1N3JsgFNBYtwc7U5pXYR9vD/B00o8JxnAnVU9/zijcRabBqKf3pAr5SeHuhbN3NSLZLbE8Id3NSFd7mDts7+LjTpIyI1mcJd5rMlhMH+7pG8Ah1ZGcHYW7yBTxqLoxapLCD/b0RnM9oYxb1Fa7hEYNuv3HrjEzKJz6ahcrjJAlD3r6kDP5SEKqzct0V+zZ0Kdd5CKZ1NM32qTbdSjHT/lDVoJFRAWjnj4i6QZ6S4EinfBjgkGUqac/Kwp3kQvosZ6+EVGbDRPrnHjgD/f03YJoUBCG8cTMdoP+gnr6s6JwF7moDvT0URcGC5On1gMg2bnf0+eEXkGyPaB295All5ERX2qopz8jCneRC84jyGc45PQPiAbgs0beMhrrMcmuEe8N8P6A8ps3YNwLsUlMvXOVZKdN1GvQX4oYzlU9fVGfvM+8GxPPlZWjU7iLyBPvQH20p09aRpm0qd9NiOP34u/chgN728BoA7Nv3iBstml1L5N0GgyWa3RXAvaEzb/V05+cwl1EjuR+T//QXbHNJulsSqOZYnc2KVfXCM0G1mgQ4hhLU8rOHv7VvyftXCXenSX0mhQzoZqhj5HXTT39KVC4i8jR3V9tk1az+N5KIJtJKBttajMpydwM5TursLlFNNvG4pioPYNnGcU7twmb27S6lylma/tnxT4in0kYLCZ0l9XTn4TCXUSeWtXT319tA/n1CKxO0YpphohydY1ifYNocaHaoTJJiBYX8G6X8qtvEi3MQTQ+rePZGZLtGaJetXtltZ7eyGZRTfMUFO4icjz3V9tERrIHWQOsiPFghLlZivUNPM+xpOpVzAxrtbBajXKvB/548e55jt9ZJ/7Qt9DMSyy0waqjCK1Eh30/BYW7iJyMjfn9kBm2xTHRbHvs35V7XcpOpwr++zdU6eTlY1GLJSIyhTRzF5GTecqZtZcl3u+Pvc6Hw9MZkyjcReT4LIO46w9OhrK8xAqn3OtWD4geLsnLfp+yswfNBhaN6W4aNezSEsVMje61GbKWMZw3HCYunZTxFO4i8vRKiHoQZVXAN9dK0q0htTtdytU1fGOLMD9HSFNgdEPTbgcPhn3Le/F2ffI692bCYKVGdyWiaBhl5NVOlyqRn4rCXUSOziFkEPWqjcTSbad+r6Bxp0e0vkv5zi0siomWFrEQcHe836fc28NWluDKEv2rLfJWPPFF16xu9JcMt+pc17JmWgJ5DAp3ETmaEuIuhBzCsJqt19YHpHe7cOsu5fYOYXb2wdJHzwvK3R08SbAPvp98uUX/Up3eciBvHT4NL+PqzFbtMXN8CncROZxXYR71HHOobzrpVknj1h5RZ0D5xjchywEoN7f2rwsBu3oJLi/QvdpiuBDTWzI8jA7xmJDvHkYHfGi2fiIKd5GLbtSfh3zyQ8yrnSEbayXpbkG0mxENCiyKCN/2gfE3JAUjW2yM6c8n7wgpp0fhLnJRPdKfx12fuFtj3KseV9/MCYVDHFXBPR9jxYSnD5DNRtVsXf35M6dwF7mIxvTnyW4OE87eCFlJyKstBBwo6oHhTMRw3iiS8WldJuCxqT8/Iwp3kYtkXH++UdBY7RJtdccfumGG1Wt4ZAwXGvSXU4qa0Vuu+nMPE6biBmVa/dJs/dlTuItMk7KqWibdNRry0a/BaG36Tk7tVgdz8F4f+oPHLzLD04RssUl/OaG/HMjVnz/3FO4i0+Bgf37YrZzuNO45Udepb+WE3CnmGlg/w68vjb/GjMFibX+1i/rzc0HhLnLeldULniEb9efrRTV7H8PdSXolUa/EqCb4RSsmu1Sr6pNx1xgM5gNF3SgjyJuoPz8HFO4i59Uj/Xlts6S2UdK42yUMx78yaoXjSUTRTCkTo7+YUNaguxxgUncOOD46gQnN1s8JhbvI8+qw3Rbvz9bv9+d3C2rrQ9K1Pbi9hg/GdOeAXVrGWk3K2YT+Ykp/KZA3R/15/ZChRKa9Xc4ZhbvI88ZHuy32mLjuvHqcU99wahs5jdUutraN31zFajVCtP+tXXarU4/s2mXKxTa9q0315xeAwl3keXKgP7fMSTsT0t0h7VR7uyRre/iNVbzXJ5qbw+Lq29qzjGJnF2s14IUrZCtNBsuju0XVn089hbvI82Bcf75ZknQm7AngTrI9xO5s4Kt3Cc0mtrBQ3WTkTrm3hw+H2AtXKFfm6F1pMlyI6C+OZuvqz6feicPdzCLgFeCmu//oyYckcsEUY/rze0PSu3vYXu/xx4/OFy03t/BOF6IIHwzxQXWKkRcFNjuDvf8Fhpfuz9YDZTq6W1R7o18IpzFz/wTwGjB7Cs8lcnE8erfomP7cH51a5zlEEba8gLXbMPf4t11o1ikWW/QuNxksBgbzAQ+j2XqCZusXxInC3cxeAH4E+E/AJ09lRCIXQTHa26WodmRsrhVH68/n2nDtCvlCvVrBMu6pGzHD+YTeSqBMjCKpgl2z9YvlpDP33wR+CWhPeoCZvQS8BBDPLZzww4mcc15VL1G/2omxcc9JN3Mat7vY3U389h1CY0x/Pqj6c1+ZpXelRdYOE/d0KVIYzlb7vuRN8OQZ/zfKc+HY4W5mPwrcdfcvmtn3TXqcu78MvAxQv3b9Kc9JF5kelkPUdUJpxF2nsVZSW+uT3Ongt1bxTpfQblcd+jDDvaTs7I3682sMV1qjvdGr/vwwRVotcdRs/eI6ycz9I8CPmdkPA3Vg1sz+0N1/8nSGJnKOOBO3ywWIhhAG1Wy9uVYQd0tab3UIg5zi5ipkGaQJ5aB/4KIIe+81yuU2vStj+vNJSxgNhbocP9zd/ZeBXwYYzdz/vYJdLqQD/flhkl2vttjdLQj9Ao8iylASPnB97OM9GIOVBsOltJqtqz+Xp6B17iLH9Uh/Xt/yiacShWz0mI0MKwEz8rkag/e2CBNm/G7Qn7f9/rwBPmFzL5FHnUq4u/tfAX91Gs8lch6M68/TjQFWjH9ZyXLHQoA4woG8GcjaEb3lQHFYYJupP5dj0cxd5Gk4RP0D/fl68WC1C+tbkI2ZuhvY/CxEEcOlBoOlGsNZYzBf3S162IujZaLVLnI8CneRg4oqvCdt2GVFdRjGfn+ek6x1CVmJt2dgrzvuKjwODC41GS7v9+dlMjpbVDNyeRco3EXg4f7cqXrxMayAxkaB5fv9edmuU3hJWYtwZsZel88k9BcCwzn15/JsKNzlwnu4P4fGekGYsF8XOPFeWb1ASrUCMm9FZDMpWWv8trkOZDOGx+rP5dlRuMvFNaY/r23m1G9XNcvEywyKmRpFK2EwGz/oz7FDTjK6P1tXfy7PiMJdplfJxKWJ5veXMFb9eWOjpHanR7y+h79zm6I7ZjdGwN5zGVteIGtFDJfix/vzw24c1YZd8gwp3GX6HOzPD0lUK6CxVlDbzKjd7mKr6/jaBmFmBltuVg8qnbLTwQ3s+nsolmfoX2mqP5fnnsJdpsqj/Xlte/Jto1Hfqd/pEd3r4G/fqt63tFitRwfK/oBydxdbWoCrywyutBgsJnRXgvpzee4p3GU6TFh/nm4MJs7dQzerDpO+t0WYbRNqteqpypJyZxcPhn3re6vZ+qUGvaVANhtwc+22KM89hbuce5ZD3B315x2nfq+kfrdHvNaFe5tU61UOGL3pu3swrE4vKrd3Htr3yy4tw9Ul+ldaDBbiam/0qDqarmjoMGl5/inc5fzy6qCLMHTCg/48p3Z7D1u9h6/dw2q1/SAuHe/3oZYSrlyq6pbo8U7FzcgXG/Qv1ekvBbKZgAcnbxqu7xg5J/SlKueSZfdn65Du7s/Wo/UO/s5tKP3h/nwwoNzZxVaW8KvLDJYblOn4PXPLpKpfusvV9rplDYq6Zutyvijc5Xwpq249GlZ1THOtpLaRka7uYbfX8XubhHabUD/Qn+/u4vBQf95fCNV5ouM+RGLVYdLBKTRbl3NKX7ZybtgQ4t5otr7t1DdLGqtdwt2d6jDpwRBrNPA8o+hk4OD9PrayBFeW6F9tMlhIHvTnE286suoko7Km2bqcXwp3eT44hIzHXvu8L+TV34esmq0nOwX1mx2sdLzXx5oNaDYevsgMu36FfLlF/3KD/mIYnT2q/lymn7685cw96M/9kGmyO7UtJ911ats5IXOKdo3Qyyivr4yfYJsxXKoxWIjVn8uFo3CXs/Nof75eVLP3MayEqFedPXp/w66yGdO7XJ84A3eDwXwgb6o/l4tHX+pyJsb15/U7XcJwwoZdhUMcKJspHhm95YQyht5KwKNDNuzCKerqz+XiUbjLsd2vx8tRnVK6TezMHygh7j3cn9fuDand7eHrm9WNRWOElUVoNshbCYPFlMGcPejPi/rkD+uRwfgVjyJTTeEuxzYYdRxbWROjpDOske2kmBvpPZtQsVSHYdS2nHTHSXdy4k4B7oTFeXx+duzH8lpM/0qTwaL6c5GjULjLsfVH4X43m+Xm7hyURrmTUluH+j3Hxs2nHdKOk+zu9+dlLSZvxWTtGPPxc/AyMnorkfpzkSPSt4ccW98TSjf6ZczdvTbZRo14x5i5Ac3bwwldiRMNRlvxWrVlbm85wWPoz4WJOywW6eiQDPXnIkeicJfDObgbRVmlbuZVgb1b1uiUKRtZi51Bg343ha2Y1g1ofbNHbXUXxs3CowiLIoo0IltqMJyN9/tzczxMSO0AeR315yJHpHCXiUJulKmT7yV0Q0k/j7k9nOeFdIMbzAPw1mCJmzuzFOs1GjeN9jeHJK/fplxbf/jJihKCYZcv4VcW6V9pMNT6c5F3jcJdJgqZkWwHvNugW0v54r06X59f4dW5q1xvbbAxmOGNO5fovdlm7uvG7Gs7RK+9QxRibHkZAM/zam+XRgN74QrZSovBSp3eSiBvBvXnIu8SfUvJoQzDMghZjHcidlfr7KRzvNG+Bl2jdiuw8npO49XbhHvbhEYTM8Pd8W6PstfFrl6Gywt0rzQZLsT0lqquXf25yLtH4S5PVkIYQjRgtEVAjL1VUtt22l/dJrz2NlEUYUkKeU7pTrnbwZoN7EMfqGbryzW6KxFFwygjyJuoPxd5Fync5YniPQhFdaB0ulNS28yprw0It+7BzTUoy4dOMSIK2LXL+MocvdHa9P6iVatdGtVpRpqti7y7FO7yREUDrAtglInjUcCDEWba2Ifajz3eDQbLjdFsPVDWjDKGvIFm6yLPiMJdnshjyNrVJl+OkTUjikaDZKUxdrmjA4N5YzAfNFsXOSMKdzma+yGdGFHX6V4aTcHLMWvZDTCjTEaz9Qk3JonIu0fhLk/FY8jbhuXVi6yTpuNlAp4+06GJyAEKd3l6Bp5AMeEMUhE5e/qBWURkCh073M3supn9pZm9ZmZfNrNPnObARETk+E5Sy+TAL7r7l8ysDXzRzP7C3b9ySmMTEZFjOvbM3d1vu/uXRn/eBV4Drp3WwERE5PhOpXM3s/cB3wl8YczfvWRmr5jZK8Xe+CPURETkdJ043M1sBvgj4OfdfefRv3f3l939RXd/MWq1TvrhRETkCE4U7maWUAX7p939c6czJBEROamTrJYx4HeA19z9N05vSCIiclInmbl/BPgp4PvN7G9Hv374lMYlIiIncOylkO7+f9BWUCIizyXdoSoiMoUU7iIiU0jhLiIyhRTuIiJTSOEuIjKFFO4iIlNI4S4iMoUU7iIiU0jhLiIyhRTuIiJTSOEuIjKFFO4iIlNI4S4iMoUU7iIiU0jhLiIyhRTuIiJTSOEuIjKFFO4iIlNI4S4iMoUU7iIiU0jhLiIyhRTuIiJTSOEuIjKFFO4iIlNI4S4iMoUU7iIiU0jhLiIyhRTuIiJTSOEuIjKFFO4iIlNI4S4iMoUU7iIiU0jhLiIyhRTuIiJTSOEuIjKFThTuZvZDZva6mb1hZp86rUGJiMj0bcklAAADxElEQVTJHDvczSwC/jPwz4BvB37CzL79tAYmIiLHd5KZ+z8G3nD3N919CPxX4GOnMywRETmJ+ATXXgPeOfD2DeCfPPogM3sJeGn05uBrv/LJV0/wMafJMrB+1oN4TuhzsU+fi336XOz70NNecJJwtzHv88fe4f4y8DKAmb3i7i+e4GNODX0u9ulzsU+fi336XOwzs1ee9pqT1DI3gOsH3n4BuHWC5xMRkVNyknD/f8AHzez9ZpYCPw78yekMS0RETuLYtYy752b274D/AUTA77r7l59w2cvH/XhTSJ+Lffpc7NPnYp8+F/ue+nNh7o/V5CIics7pDlURkSmkcBcRmULPJNy1TUHFzK6b2V+a2Wtm9mUz+8RZj+msmVlkZn9jZn961mM5S2Y2b2afNbOvjr4+/ulZj+msmNkvjL4/XjWzz5hZ/azH9KyY2e+a2V0ze/XA+xbN7C/M7Ouj3xeO8lzverhrm4KH5MAvuvs/BL4b+LcX+HNx3yeA1856EM+B3wL+zN3/AfAdXNDPiZldA34OeNHdP0y1WOPHz3ZUz9TvAz/0yPs+BXze3T8IfH709hM9i5m7tikYcffb7v6l0Z93qb6Br53tqM6Omb0A/Ajw22c9lrNkZrPA9wK/A+DuQ3ffOttRnakYaJhZDDS5QPfPuPv/BjYeeffHgD8Y/fkPgH9+lOd6FuE+bpuCCxto95nZ+4DvBL5wtiM5U78J/BJQnvVAztgHgDXg90YV1W+bWeusB3UW3P0m8OvA28BtYNvd//xsR3XmLrv7bagmiMClo1z0LML9SNsUXCRmNgP8EfDz7r5z1uM5C2b2o8Bdd//iWY/lORAD3wX8F3f/TmCPI/7oPW1GffLHgPcD7wFaZvaTZzuq8+lZhLu2KTjAzBKqYP+0u3/urMdzhj4C/JiZfZOqqvt+M/vDsx3SmbkB3HD3+z/FfZYq7C+iHwC+4e5r7p4BnwO+54zHdNbumNlVgNHvd49y0bMId21TMGJmRtWrvubuv3HW4zlL7v7L7v6Cu7+P6mvif7r7hZyhufsq8I6Z3d/576PAV85wSGfpbeC7zaw5+n75KBf0xeUD/gT4+OjPHwf++CgXnWRXyCM55jYF0+ojwE8Bf2dmfzt6339w9/9+hmOS58PPAp8eTYDeBH76jMdzJtz9C2b2WeBLVKvL/oYLtA2BmX0G+D5g2cxuAL8K/Brw38zsZ6j+8fuXR3oubT8gIjJ9dIeqiMgUUriLiEwhhbuIyBRSuIuITCGFu4jIFFK4i4hMIYW7iMgU+v+P1k6OKyeCXgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFqBJREFUeJzt3VuMXWd5h/H/ix07BJyEgyHUk9jBijJFUWnSqAWiItQQiZaIcEGlIIGiCsmq1HKWUOCC3FVcIEQuKqRpwqFKlKo1kUAU0VihqKpaRc0BlQSbQ8AOQ2xi1BIs1NqkvL3Ye3n2LO/DOnxrfd/61vOTIs+MZ89e3oqfef3utdeYuwsAkJcXxT4AAEB4xB0AMkTcASBDxB0AMkTcASBDxB0AMrQy7mb2eTN7zsyenPnYy83siJn9YPrry7o9TABAHVUm9y9KelvpY3dKetjdr5H08PR9AEAirMqLmMzsgKSvuft10/e/J+kt7n7SzF4j6Vvufm2XBwoAqG5nw9u92t1PStI08K9a9IlmdkjSIUmyi3b93q69Cz8VADDH2Wc3f+7ue+vcpmncK3P3DUkbknTxvit9/59/pOu7BICsfP+THzlR9zZNz5b52XQdo+mvzzX8OgCADjSN+1cl3TF9+w5JXwlzOACAEKqcCvmApH+XdK2ZbZrZ+yR9StItZvYDSbdM3wcAJGLlzt3d373gt24OfCwAgEB4hSoAZIi4A0CGiDsAZIi4A0CGiDsAZIi4A0CGiDsAZIi4A0CGOr9w2Cpnrzq37f3dz+yKdCQAkI9ocS9HfdnHCT4A1BNlLbMo7Ms+v+5tAGDMeo277/JWkSbyAFBN9J17E+zpAWC5Qca9bDb2hB4AMon7LEIPAJmf586OHsBYZTe5z8M0D2BsRhH3WYQewBhkvZZZhbUNgFyNbnKfh2keQG5GPbnPwyQPIAdM7nMwyQMYOuK+AqEHMESsZWrgCVgAQ0HcGyDyAFLHWqYFVjb17TnhnX3tM/uts68NDA1xD6QIPZHf0mXI694f4cfYEPfAxh75voNe1aLjIvrIFXHvyJhWNqkGvYp5x07wkQPi3oOzV53LMvBDjvoy5T8XsccQEfee5DTJ5xr1RYg9hoi4RzDk0I8t7PMQewwBcY9sSE/AEvb5Zh8XQo9UEPdEpL6XJ+zVEHqkgrgnJNV1DWFvhvUNYmp1+QEz+7CZPWVmT5rZA2Z2cagDG7tULnFA2MPZc8LP/wd0rXHczWyfpA9IutHdr5O0Q9LtoQ4ME6lEHmERenSt7Vpmp6QXm9mvJV0i6dn2h4R5YqxsugrPZU+fnfvx5w/u7uT+UseeHl1oHHd3/6mZfVrSM5L+R9JD7v5Q+fPM7JCkQ5K04xWXN707zEj9yVdpccBD3CbnbwKEHqE0jruZvUzSbZKulvQLSf9gZu9x9/tmP8/dNyRtSNLuq9f4N2ggXU/yTab2JkFvYt795Bh8Qo822qxl3irpx+5+WpLM7EFJb5J039JbIbiY58r3FfRVcg9+EXoij6raxP0ZSW8ws0s0WcvcLOnRIEeFRvqOfCphXyTH4DPNo6o2O/dHzOywpMclvSDpCU3XL4irj8inHvZFysc95NgTeizT6mwZd79L0l2BjgWBtdnLn9lvC/fuQw37PLN/FkKPnPAK1ZFoMs3PC3yVsO86tnn+7XPra5XvL7Zcpnr285CI++iEXNnMRrzJ56Qe/qFP9Uzz40bcR6rqymbReqZK2FeZ9zVSDT6hx9AQd6wM/bzAn1tfCxL4svLXTDH2uYSeyOeNuGOb8nVsitgXITizf3JtuD0nXDp4UJc9fbaTyBdSj/2QQ880nzfijqUWXbTs7FXFW7t1mbqZ4udJOfaEHikh7pjrwNrppb9/yxXHdOTUuo7/4V6d2X+x9vQwxc+T6pk5OYSeyA8bcR+hVeEuu+WKY+ffftelj297+/AVN+jI2rqOXzWJ/FqPU3xZqlP9UEPPND9sxH0E6sZ81qKwL/L8wX7XNMukONUPPfREfjiIe6baBF3aHvVliuAfkXRceyXtUrGHnxU79ilO9UXohxh5idCnjrhnom3MZ80Le5WpXSqegJ0EXtoKWDmmKcU+duiZ5tEF4j5wXUe9jgNrp3V8c+/5wO854dtiNRuxrs6Tb4LQt8M0nybiPlB9Rr3K1F6cPVMEXtr6i1785S9iVZ7mU4m8ROjbYppPR+MfkI04DqydTmpaX2T2/Pgz+23bX/ZyqGJHdJFdxzbP/xfbZU+fHdTVOPnh3/ER94EIHXWpWtirTO3F5xRfrzjO8gugyoGfjfy59bVkIy8pucgPJfREPh7WMgMQI+qhlH+Yd/k6Nc8f3H3BLr6QQkzLWNs0w16+f8Q9YaGjLtULe9UzZOaZ3b2XLzM8L/DShddTJ/TVDem0Svby/WAtk6ihhn32Ppb9Geb9xV4WJtY21QxpbcPKplvEPUGxwx7S7J9l3g6+HPnyLr4s5cBLaT4Jm3roiXw3iHtiUgh7k6l99jbL7m/eVSZzm+ILKYY+ZUQ+LOKekBTC3oXyn6tO4Ic8xc9KLfIph57Ih0HcE9FF2FPSNPBSHlN8gWm+OiLfDnHPWJOpvc0ZMnXvu27gc5niC6mEPvVpnsg3Q9wTkOs6pmzen7NO4KW8pvhZqYU+RQS+HuIeWUrrmLZTe/n2bb7BNA28NOzIS2ns51Od5pniqyPuGUpxai9Und6l+adKFlataaR8Ip9K6FNC5Fcj7hGlNLV3peo3mkWBX6XKKzKHHnkpjdCnOM0T+MWIe2ZSntoLi76pLZvgl6n6kvscIi+lFfoUMMXPR9wjSW1qD3WWTFVdBH5skZfi7+dTizy2EPcIugp7qlN7qONatoMv1LlwVo6RjxX6VCLPFL+FuCOoOv8CqDu9F0JO8VJekZfihj6lyI8dcUeS2gZeqn/529wiL8Vb26Tw5OvYA0/ce5biSqaPffui41v2eMQIvDTMV7uuMtZpfsxrGuKO6PoIPFP8lrFGfmxaxd3MLjezw2Z2zMyOmtkbQx1YjlI7Q6ZvXT3hW/Un+jSd4nOPfN+hjxX5sQW+7eR+t6RvuPu6pNdLOtr+kFBXamfJNFnztJnepW4DL+UdeSnONB8j8mMKfOO4m9mlkt4s6V5Jcvdz7v6LUAcGzEoh8FKe+/hZMaZ5At+NNpP7ayWdlvQFM3vCzO4xs5eUP8nMDpnZo2b26P+d+VWLu0Pu+lxbNdnDF3Kf4gt9Rr7vKX4MgW8T952SbpD0OXe/XtKvJN1Z/iR333D3G939xh17Lmj/aIx9315os0IKOb0X2k7xRD6sPiOfe+DbxH1T0qa7PzJ9/7AmsUeP2u7b+77swCqrvgmmFngp/1VNoc+VDYFvr3Hc3f2UpJ+Y2bXTD90s6btBjgpoicB3q4/I9zXF5xr4tmfLvF/S/Wb2n5J+V9JftT8k5K7tvzaqXh6YwHevr8ijvlZxd/dvT/fpv+Pu73T3/w51YDkZ4769zbqnyuOVcuDHHPmuQt/1FJ/j9M4rVIGStoGXxjnFF7qOfFdyCzxxRxSrVjMxp3eJwIfQVeQJfDXEfcBSe2VqypoEPoSxB17qJvLs4Vcj7hi0pj97tYoQ07tE4Auh9/JdBT6X6Z24I1mhn4iOtZ6RCHxZqMgzwS9G3IEVCHx3Ug18DtM7cUc0oZ4zqLOaibV7x2Ihpngm+AsR946N8Rz3HDG9d69t5An8dsQdnQlx3ZouvjnGnt4J/HKxflJU2dBXM8QdWejyrJlCqOkd1TQJPNP7FuIORMD0Xg2Bb464A0hazBXNkFczxB2jFHvvDnSNuAM1hNy7s5qpru70zmqGuAMYiBTOoBkS4g4AGSLuAJAh4g4AGSLuAAaBJ6DrIe5ADSHPwuAJwu7wamLijpEa8otTxoipvT7iDiBpMcM+5Be7EXcgAlYy1TQJOyuZCeKOLOx+Zlfn98GrHvvFKqadnbEPAPk6/MsbYh/CXLH37Uzty7WJesipfcgrGYnJvXPHN/fGPoRBS+XxCzW1E/blUgl7DpjcEc2RU+tBvk6dlUzsqR3ztV3BEPYLEXdgBab27oTYq3cR9qGvZCTijoSFXsk0mdoJezdCPVnKxL4Ycceg9XGWTFuEfUvIM2C6CnsOU7vEE6qDFmpnPQaxpnbCPnFufW0QYc8JkzuiWPWNqcpKpurUTtjj6Oo89S7DnsvULhF34AKEvZ0uX3xE2KtjLdODVM7V7lPXL2Dqamon7M2FXr2UsYqph8kdvWu7kkk17GOMeh+XCOgj6rlN7VKAuJvZDkmPSvqpu9/a/pCA9gh7t/q67gthby7EWuaDko4G+DqIILXrv4Sa2usg7NV1vXopPH9wN2FvqVXczWxN0tsl3RPmcFDXmE6H7GId0ybsu45tjiLsRdBzmtalvMMutV/LfFbSxyTtWfQJZnZI0iFJ2vGKy1ve3XAd39yrA2unYx9GdMu+GYV44rnPsOeu70vu9vmEae5hl1rE3cxulfScuz9mZm9Z9HnuviFpQ5J2X73GVZvQSJWpvWrYifpyOUddGkfYpXaT+02S3mFmfyLpYkmXmtl97v6eMIeGoWqyx182tacQ9tyjHuMHY8Q4tXEsYZda7Nzd/ePuvubuByTdLumbhH25rs53H8revavjJOzN9blLL/T1ZGnZmMIucZ47EtBmaq8SdqK+XawfXxfzRUhjC7sU6BWq7v4tznGPq81U3MfpkE2OL1bYcz0LJsaULsWb1AtjDLvE5N47zprZbtHUHiPsuQY9ltiXCxhr1AvEHdE0fQ4idNiJelixoy4Rdom4R9HV9H7k1LpuueJY8K9bx7wVT92VzLKpfVXYxxz1sQe9QNgniDskTaL8rksf7+3+mqxjQoU9p6jHDLpE1FNG3CPJeXpvqmnYxxb12EGX0oq6RNjnIe7o1LyVzLypvcuwE/VwiPpwEPeIxji9130StU3Yc4g6QZ+PqK9G3CNL6dTItnv38pOpVZ9InTe1jznqKQRdSjPqEmGvirgnoIvApzi9V13HjDHsBH01ol4Pcc9YzMBXmdpDhn2IUU8l6BJRzxFxT0Qq65muToksT+11wk7Uu5Fy0CWi3hZxT0iu6xnCnk7QJaI+FsQ9MSkEvsn0Pvtk6rKVTNWw5xB1gl4PUQ+LuCcohcCHMju1Vwl7Dk+YphL1IQRdIupdIe6Jih34prv32al92TntuU3rqQRdIuqYIO4JK+IYMvJ9BF5avmevGvbUgy6lE/WhBF0i6n0h7gMQeorvc0VTXsWsWsMQ9OqGFHSJqPeNuA9E6Cm+WJ+sinyV6b14MrX4msWxLgv70Cb1VIIuDSvqBD0e4j4wQ5ni64Q91agT9OaIenzEfYC6mOKXBb7O7r28ax9a1Al6O0Q9HcR9wGZD2jb0Vdc0q24vTab2IuqLzoBJKewpBV0i6giDuGciVOjbRP745t7zYU99p07Q2yPoaSPuGSqvRprEvhz5RauZ8mV+i7CnFPICQQ+DqA8DcR+BkOubeY6cWj8/tUueVNgJehgEfXiI+8is+klI5fjP7tLnrWyOnFrXyX/bp1eecF329P9GD3tqMZeGG3SJqA8Zccc2VX4M3t+Urhez54Rr7z8+3eVhLZRizAtDjTpBzwNxx0qLfnj1nkjTOkHvBlHPC3HHNotCLumC0xv7inrKMZcIOtJE3LE06NL8FyJ1GfbUYy4RdKSPuI/Yqqgvc259LUjghxByadgxLxD1cSHuI1Qn6st+YHU5zMtiP5SIzyLoGDLiPhJNpvR5YX/+4O6FlxQYYsDLCDpyQdwz1mbtsmxiXxb4ockh5hJBx4WIe4baRL2qoQY+l5gXiDoWaRx3M7tS0t9KukLSbyRtuPvdoQ4MzfQR9iHJLeYSQUc1bSb3FyR91N0fN7M9kh4zsyPu/t1Ax4aKYgW9HM7Yk3yOIS8QdNTVOO7uflLSyenbZ8zsqKR9koh7T7qM+pn9tnTvPs9sXLsOfc4hLxB0tBFk525mByRdL+mROb93SNIhSdrxistD3N3oDWH1sii+VaI/hnAvQtARSuu4m9lLJX1Z0ofc/Zfl33f3DUkbkrT76rV6oyC26TvqTab3VcYc7kUIOrrQKu5mdpEmYb/f3R8Mc0goG8KkjnoIOrrW5mwZk3SvpKPu/plwh4RCClHvYnofK4KOPrWZ3G+S9F5J3zGzb08/9gl3/3r7wxqvFIJeRuCbI+iIpc3ZMv8qif9zA0kx6rMIfHUEHSngFaqRpR71WQR+PmKOFBH3SIYU9VkEnphjGIh7j4Ya9LIibmOJPDHHEBH3HuQS9bKcI0/QMXTEvUO5Rr1sNoRDDD0hR46Ie2BjCfoiKU/zRBxjQtwDGXvUyxaFtOvoE3Bggri3QNDrI75AP4h7A0QdQOqIew1EHcBQEPcKiDqAoSHuCxB0AENG3EuIOoAcEPcpog4gJ6OOO0EHkKtRxp2oA8jdKOJOzAGMTdZxJ+oAxirLuBN1AGOXRdyJOQBsN8i4E3MAWK7XuNu5dlcEJOoAUE3vk3sR6LNXnat9GwBANdHWMssiT8wBoJ3oO3dCDgDhvSj2AQAAwiPuAJAh4g4AGSLuAJAh4g4AGSLuAJAh4g4AGSLuAJAh4g4AGSLuAJChVnE3s7eZ2ffM7IdmdmeogwIAtNM47ma2Q9JfS/pjSa+T9G4ze12oAwMANNdmcv99ST909x+5+zlJfyfptjCHBQBoo81VIfdJ+snM+5uS/qD8SWZ2SNKh6btnv//JjzzZ4j5z8kpJP499EIngsdjCY7GFx2LLtXVv0Cbu836skl/wAfcNSRuSZGaPuvuNLe4zGzwWW3gstvBYbOGx2GJmj9a9TZu1zKakK2feX5P0bIuvBwAIpE3c/0PSNWZ2tZntknS7pK+GOSwAQBuN1zLu/oKZ/aWkf5K0Q9Ln3f2pFTfbaHp/GeKx2MJjsYXHYguPxZbaj4W5X7AmBwAMHK9QBYAMEXcAyFAvcecyBRNmdqWZ/bOZHTWzp8zsg7GPKTYz22FmT5jZ12IfS0xmdrmZHTazY9P/P94Y+5hiMbMPT/9+PGlmD5jZxbGPqS9m9nkze87Mnpz52MvN7IiZ/WD668uqfK3O485lCrZ5QdJH3f23Jb1B0l+M+LEofFDS0dgHkYC7JX3D3dclvV4jfUzMbJ+kD0i60d2v0+RkjdvjHlWvvijpbaWP3SnpYXe/RtLD0/dX6mNy5zIFU+5+0t0fn759RpO/wPviHlU8ZrYm6e2S7ol9LDGZ2aWS3izpXkly93Pu/ou4RxXVTkkvNrOdki7RiF4/4+7/Ium/Sh++TdKXpm9/SdI7q3ytPuI+7zIFow1awcwOSLpe0iNxjySqz0r6mKTfxD6QyF4r6bSkL0xXVPeY2UtiH1QM7v5TSZ+W9Iykk5Ked/eH4h5VdK9295PSZECU9KoqN+oj7pUuUzAmZvZSSV+W9CF3/2Xs44nBzG6V9Jy7Pxb7WBKwU9INkj7n7tdL+pUq/tM7N9N98m2Srpb0W5JeYmbviXtUw9RH3LlMwQwzu0iTsN/v7g/GPp6IbpL0DjM7rsmq7o/M7L64hxTNpqRNdy/+FXdYk9iP0Vsl/djdT7v7ryU9KOlNkY8ptp+Z2Wskafrrc1Vu1EfcuUzBlJmZJnvVo+7+mdjHE5O7f9zd19z9gCb/T3zT3Uc5obn7KUk/MbPiyn83S/puxEOK6RlJbzCzS6Z/X27WSJ9cnvFVSXdM375D0leq3KjNVSEraXiZglzdJOm9kr5jZt+efuwT7v71iMeENLxf0v3TAehHkv4s8vFE4e6PmNlhSY9rcnbZExrRZQjM7AFJb5H0SjPblHSXpE9J+nsze58m3/z+tNLX4vIDAJAfXqEKABki7gCQIeIOABki7gCQIeIOABki7gCQIeIOABn6f8Bln6LeuNr0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0512718324071463\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAF41JREFUeJzt3WuI3Xedx/HP18REW9PUS7SaaZtapLMi67ZbdtWyItaCi2J94EIFpYhLWNj1DlJ9oM8WH4jog0XI1ttiUSQWFFd2DXVlWXYpG1tZq4mXaBJHExvZtQ1FE6vffTDnnzlzci7/++/2fkHpzGTOzG/+ybzPd37nf/5j7i4AQF6eEnoBAID+EXcAyBBxB4AMEXcAyBBxB4AMEXcAyNDKuJvZp83sUTN7ZOptzzKzI2b2o8n/nznsMgEATdSZ3D8r6bUzb7tH0gPu/iJJD0xeBwBEwuo8icnMDkj6mru/ZPL6DyS9yt3PmNnzJX3L3W8acqEAgPp2trzd89z9jCRNAv/cRe9oZgclHZQke+quP921b/u7+q50niFrFy30EgAU6MIvNn7l7vua3KZt3Gtz90OSDknS0/Zf69f/zXsv/dmF6y4O/el7t/v0rtBLAFCYH37ovaea3qZt3H9pZs+f2pZ5tMmNU4w6AKSk7amQX5V09+TluyV9pe4NCTsADK/OqZBfkPRfkm4ysw0ze7ukj0i6w8x+JOmOyesrpbS/DgApW7kt4+5vXvBHt/e8liRcuO4i++4AosczVAEgQ8QdADJE3FvgQWEAsSPuAJAh4t4S0zuAmBF3AMgQce+A6R1ArIg7AGSIuHfE9A4gRsS9BwQeQGyIe08IPICYEHcAyBBx7xHTO4BYEPeeEXgAMSDuAyDwAEIj7gMh8ABCGvwXZJeMX+zRvz2n0vhtXuevt9BLQOGI+8AIfDOpxHuVpl8HdwboG3EfAYGfL5eQ96HOseAOAE0Q95EQeGLe1bLjR/gxi7iPqMTAE/RxLDrORL9cxH1kpQSeqMeB6JeLuAeQc+CJehrm/T0R/LwQ90ByCjxBz8Ps3yOxTxtxDyj1wBP1vBH7tBH3wFINPGEvD7FPC3GPQHWpglQiT9ghbf93QOjjQ9wjksIUT9gxD6GPD3GPTMyBJ+yog9DHgbhHKMZtGsKONgh9OMQ9YrFM8YQdfSD04+J67pELfV14wo4h7Dnl/NsaGHFPQOjAA0Mh8sNhWyYRIfbhQ3zT7T1xYfTP+diNu0f/nNiOLZv+EffExLIP34cQIZ+nzjq4AxjPnlNO4HvQKe5m9h5Jfy3JJX1X0tvc/bd9LAyLjTHF9z21xxLytlatn/j3q/r3R+Tbax13M9sv6Z2SXuzuvzGzL0m6S9Jne1obVoh5ik895k0t+nqJfjdEvr2u2zI7JT3dzH4n6QpJv+i+JDQR2znxpUV9lXnHg+A3R+Sbax13d/+5mX1U0mlJv5H0DXf/xuz7mdlBSQclacezr2776bBCDFM8Ya+H4LfHfnx9rU+FNLNnSrpT0g2SXiDpSjN7y+z7ufshd7/V3W/dsefK9ivFSheuuxjktMm9Jy4Q9o6qY8ixXI3TJ+vpsi3zGkk/dfdzkmRm90t6haTP97EwtDfWVg0RGs7ssWWyvxxT/HJd4n5a0svM7AptbsvcLuloL6tCL4bcqiHs4yL287EXv1iXPfcHzeywpIckPSnpYUmH+loY+jHEFE/YwyP22zHFX67T2TLu/mFJH+5pLRhQ08ifv97Y10zIdOxLDT1T/HY8Q7UwTSI/L/DLpvZdxze6LU7SxfW1zh+jdKWHnil+E3EvVN3IL5vg+4h5049J/JspNfQEnrgXr+2e/BBh7/p5Cf9ypYW+9MATd0haHvnp6f2xG3dr74kLuri+Fizwi8xbD8Gfr5TQlxx44o5tZp8EVcV++hvk/PVP24z9jTdq74kL0UV+GsFfLffQlxp44o6lFj3j9cJ11Uu7tVfxTfHLEPzFqtDnFvkSA0/cccmBtXO13/eOa47ryNl1nfyLfZuTfAJT/DKz6y499jlO86UFnrhnpkmg67rjmuOXXn7TVQ9te/nwNbfoyNq6Tl63Gfm1xKb4RYj9lpym+ZICT9wTM0S8l1kU9kUeuzG9bZo6iH0+kS8l8MQ9YmOHfNp01KXFYa/efkTSSe2TtEvVPvw0Yp+PHCJfQuCJe2RCBl26POpNbT4Auxl4aSsEs/HLOfalhD71yOceeOIeWOiYS6uDXmc7Rtr8Wk5u7LsU+D2nfNs3/vSDdDGeJ9+X0qb61COfK+IeSApRb/qxjpxdvxR4aevc+OknQEmXT/O5Rr5SylSfYuRznt6J+8hSi3rdqX3W9LXk50W+lCl+VgmhTy3yuQa+9a/ZQ30H1s5d+i+0Pqf1SnUHUH3s6uucfQLU9DfQYzfu3vbNf3F9LdvYLbLr+Mal/3KU0q8MzPHy1kzuA4oh5pU2UW87tU+b/W1Qs1eZnDfFV3KN3jw5T/R7T1xIZorPCXEfSCxhbzupdw379N777EXJ5gVeuvxa8YQ+n9CnsFWT2/YM2zI9i2X7RQoT9unPuew4zPsmWvaNn0vkmspt6yb2rZqctmeY3HuUetSHNm+LRlKtKV4q5+yaRXKa6NmqGR6Tew9ymNYrbaf26dstm97nXWWyzRSfety6ymGij3WCz2V6J+4dxRT1UGFfpUvg2apZLeXIx7pNk0PgiXsHMYU9NXUDL7EXX1fK03yMgU8dcW8pt7D3PbXPrmve8eoz8ER+uxRDH9sUn/r0TtxbyC3sY6l73JYFnsg3l2Lk0R1xbyjHsA+1115njYt+jd+y841XnWVB5OdLKfKxBD7l6Z24N5Bj2PtU506i7vaM1C3wEvvxi6SyZRNL4FNF3BPTd9iHmtqbaht4pvhuYo98DPvwqU7vxL2mGKb2WCf2ZeatuemxXPWU8LpTPJFfLIXIoxniXkOuYY9laq8smt6lfgIvEflVYo58yMCnOL0T9xVyDfuYmkzvXQNP5PsRa+SZ4Osj7oUaampv8nGHCLzU7MqDRH65GCMfKvCpTe/EfQmm9nDGDLzEmTWrxBZ5JvjViHvEhgp7iL32RV9L2zvQoQJP5JeLKfIEfjnivkAMU3vplk3v0jCBl4h8HbFEfuzAp7Q10ynuZna1mR02s+NmdszMXt7XwkpX0nbMsjvSvgJP5IcRQ+SZ4OfrOrl/QtK/uPu6pJdKOtZ9SeHlPLXHdvrjmNr+cggCv1royI8Z+FSm99ZxN7OrJL1S0qckyd0vuvuv+1pYyVKf2hfdgbT9uvqY3isEflihI48tXSb3F0o6J+kzZvawmd1rZlfOvpOZHTSzo2Z29Pfnn+jw6cYRempPPextdT3uYwWeyNcTIvJsz2zXJe47Jd0i6ZPufrOkJyTdM/tO7n7I3W9191t37Lms/RhRylsyq6Z3aZzAS0zxTeQa+BS2ZrrEfUPShrs/OHn9sDZjDyy07CeTVdM7gU/T2FM8E/ym1nF397OSfmZmN03edLuk7/eyqkDYkilP18AT+frGjDyB7362zDsk3Wdm/yPpTyT9ffclAYv1Pb1L3QIvMcU3xQOu4+gUd3f/zmQ//Y/d/Y3u/n99Law0Q0/tqey39/XTU9PAd0Xgmxljih96eo99351nqE6E3pIpSdc7sjrTe1Ndp3eJwLeReuBjRtwxiC4/KYSa3gl8GJwbPwziHgEeSI0HgQ9nqMCXOr0TdySp7tbM2Hvv6Ca1CT7mfXfiXoBUHkyNBdN7WENs05Q4vRN38WBqCKu2our8ncQ+vRP4blKb4mND3IE5+pje0V2fgS9teifugfFg6jiY3tPFBN8OcQcW6Gt6J/DdEfjmiDuSNsQTmhCnPgJf0tYMcUe0eKAbKYj1dEjiDoyArZl+sD1TH3FHMdo8qMpZM/Eh8PUQdwDIEHEHkJwu03spD6oSdwDIEHEHgAwRdwDIEHEHkJwup5aWcgYUcQeADBF3AEnhCWH1EHcUo83TxPs6bY4n3vSDsNdH3BGtkxv7Qi8BWCnWX+VI3JG03ad3hV4CRtLH1F7Kg6kScQcWYksmHmzHNEfcUYRYL8uK1Qh7O8Q9sCNn10MvAXMwtcehz7CXtCUjEXcE0sedWt399lBTO2Hvhom9G+Ku/M/KOPz4LaGX0FjIv5NSrhoYs77DXtrULkk7Qy8AGBJTe1pSm9ZjPQ1SYnKPAvvuzQ11CmQfUzthb2eosJc4tUtM7hhIl62gvrZkmk7thD2M1Kb1VDC5T+S+7x6Trj+pxPrEJcLe3NBhL3Vql4h7NIbemknxQdUuxp7aCXszF9fXkg97zPvtEnFHZFb9BFVnaifscWMbZhyd425mO8zsYTP7Wh8LComtmfIQ9vGMMa1XSt6OqfTxgOq7JB2TdFUPH6toR86u645rjodexqCWbT+NPbV3CTtRr2/sSZ2wb+o0uZvZmqTXSbq3n+VgSCnvuxP2NOUa9tj326Xu2zIfl/R+SX9Y9A5mdtDMjprZ0d+ff6Ljpxte6K0Zznlvh7DHZcwtGMzXelvGzF4v6VF3/7aZvWrR+7n7IUmHJGn3DWtcmq8AbX5CWHan2uepj23DTtTrCRl0tmO26zK53ybpDWZ2UtIXJb3azD7fy6oCY3ofxlBfV92pnbAPJ/SkPmbYU9iSkTrE3d0/4O5r7n5A0l2Svunub+ltZYUbKoQx7rt3mdqHDPuu4xuEvYbQ2y9M7PNx+YEFTm7s04G1c6GXUbQ+wt426lgtdNSl8cOeytQu9fQkJnf/lru/vo+PhS05Te+Lvpa2W2CEPZzQWzAVJvbleIbqEqH33qV8999XWTa1DxF2tmBWiyXqEmGvg20Z9KrJTwWL7jzHDDtBXy2WoFdChT2lLRmJyX2lXKf3Mbdmmqy/S9j3nrhA2HsU06ReYWKvj8m9hhgeXM3t0gRN7zTrhL0uor5YbDGfFjLsqU3tEpN7Uvqe4MeY3uetuel2TF9hZ199sRin9GlM7M0R95pi2J6R8n2AdciwE/XFYo+6FD7sKU7tEnFvJMfAj31a5Lxj2CbsdfbXifpiqUQ9dNhTRtwbIvD1P06XNa4K+zJEfb4q6LFHXQo/rVdSndol4t5KjoEfQ92pvW3Yifp8qQS9EkvYU0fcW8ot8H1vz9RZV5Owr9qGIerbpTSlT4sp7ClP7RJx7yS3wA9p9lg1DfsiTOvbpRh0if31IRD3jmIKfNfID/XgatuwM63Xk+qUXokx6qlP7RJPYupFFa/QT3SSwj3ZafqOYdmdTJOwL0LUN6Ua80qMUZfyCLvE5N6rmKb4tvqe3qePSZ2wM60vl/qUXok17Dlhcu9ZDJcqkLYC32aKP/z4LXrTVQ91+rzS8js7pvX6Ug/5tNijnsvULhH3QcS2TSM1j3yXwEvL99nrhp2g54Wwj4u4Dyi2yI8deOnyrZh52zDTCHp+Yo96roj7CGKJfJetmmWqffrq41df77KwM6lvyjXoUlpRz21ql4j7qFKMfNvpvUnYS4t6zkGX0oq6lGfYJeIexPR+dMjQDzHJz+61E/X8Y15JLeq5I+6BxRD6VZGvM71PnyWz+/SuS1FfdAZM7mEvJeiVVMOe69QuEfeozE69Y8e+6yR/cmPfpbCXtqdeWswrqUZdyjvsEnGPWqjYz4v8oul99klPVdhzDnmFoKcr97BLxD0py54UNET4654+eeTs+qWpXfJsw15qzCs5RF0qI+wScc9G20sfrLpTmN5LnzfRHzm7rjP/uV/POeXae+K3WYW99JhXcom6VE7YJeJevDZ3Cv84c72YPadc+/75RJ/LCoKYb5dT1KWywi4Rd9Sw7HecpjqtE/L5cgt6pbSwS8Qdc6z6hdXVmTCpRJ2QL5dr0Cslhl0i7phYFHRp/hORYgw7EW8m96hL5YZdIu7FWxb1ZS6urwULPBFvr4SgV0oOu0Tci1U36ot+r6m0PbJ9h56A96ekoFdKD7tE3IvTZFKfF/bHbtw999mnxDguJQa9Qtg3EfdCNN1+WTaxIz4lx7xC1Lcj7plru6e+zKLpHeMi6FsI++Vax93MrpX0T5KukfQHSYfc/RN9LQzdDRH2CoEfHzGfj7DP12Vyf1LS+9z9ITPbI+nbZnbE3b/f09rQ0pBRn0bgh0XMlyPqy7WOu7ufkXRm8vJ5Mzsmab8k4h7QWGGvVAEi8t0R8/oI+2q97Lmb2QFJN0t6cM6fHZR0UJJ2PPvqPj4d5hg76rOIfDOEvB2iXl/nuJvZMyR9WdK73f3x2T9390OSDknS7hvWOAWjZ6GjPovIX46Qd0fUm+sUdzN7qjbDfp+739/PklBXbGGfNhu0EmJPxIdB2NvpcraMSfqUpGPu/rH+loQ6hg77+eut13Pdc4k9AR8PUe+my+R+m6S3SvqumX1n8rYPuvvXuy8Li8Q8rTcxL5Khg0+440DU+9HlbJn/kMTfwojGDnvf0/sqxLVsRL1fPEM1AblM68Asgj6cp4ReAJYLHXa++TCE89cb/7YGRtwjFjrsFb4J0ReiPh62ZSIUS9Snjb3/jnwQ8zCIe2RiDHuFwKMugh4ecY9IzGGvEHgsQtDjQtwjkULYKwQeEjGPHXEPLKWoTyPwZSLo6SDuaK36Rify+SLm6SLugaQ6sc9D5PNAyPNC3APIKezTiHw6CHn+iPvIcg37tOlwEPrwCHmZiPuISgj7LEI/HiKOacR9JCWGfdZsfIh9cwQcdRH3ERD2+Yj9doQbfSLuAyPs9S2LW8rhJ9oIgbgPiLD3p2sgm9w5EGPkgLgPhLDHhWCjNFzPfQCEHUBoxL1nhB1ADIh7jwg7gFgQdwDIEHHvCVM7gJgQdwDIEHHvAVM7gNgQ944IO4AYEfcOCDuAWBF3AMgQcW+JqR1AzIg7AGSIuLfA1A4gdsQdADJE3AEgQ8S9IbZkAKSAuANAhjrF3cxea2Y/MLMfm9k9K9//ojH5AsAIWsfdzHZI+gdJfynpxZLebGYvrnNbAg8Aw+oyuf+ZpB+7+0/c/aKkL0q6s+6Nd5/elVzkU1svgHJ1+QXZ+yX9bOr1DUl/PvtOZnZQ0sHJqxd++KH3PtLhc+bkOZJ+FXoRkeBYbOFYbOFYbLmp6Q26xH3er5P3y97gfkjSIUkys6PufmuHz5kNjsUWjsUWjsUWjsUWMzva9DZdtmU2JF079fqapF90+HgAgJ50ift/S3qRmd1gZrsk3SXpq/0sCwDQRettGXd/0sz+TtK/Stoh6dPu/r0VNzvU9vNliGOxhWOxhWOxhWOxpfGxMPfLtskBAInjGaoAkCHiDgAZGiXuTS9TkCszu9bM/s3MjpnZ98zsXaHXFJqZ7TCzh83sa6HXEpKZXW1mh83s+OTfx8tDrykUM3vP5PvjETP7gpk9LfSaxmJmnzazR83skam3PcvMjpjZjyb/f2adjzV43LtcpiBDT0p6n7v/kaSXSfrbgo9F5V2SjoVeRAQ+Ielf3H1d0ktV6DExs/2S3inpVnd/iTZP1rgr7KpG9VlJr5152z2SHnD3F0l6YPL6SmNM7p0uU5ATdz/j7g9NXj6vzW/g/WFXFY6ZrUl6naR7Q68lJDO7StIrJX1Kktz9orv/Ouyqgtop6elmtlPSFSro+TPu/u+S/nfmzXdK+tzk5c9JemOdjzVG3OddpqDYoFXM7ICkmyU9GHYlQX1c0vsl/SH0QgJ7oaRzkj4z2aK618yuDL2oENz955I+Kum0pDOSHnP3b4RdVXDPc/cz0uaAKOm5dW40RtxrXaagJGb2DElflvRud3889HpCMLPXS3rU3b8dei0R2CnpFkmfdPebJT2hmj9652ayn3ynpBskvUDSlWb2lrCrStMYcecyBVPM7KnaDPt97n5/6PUEdJukN5jZSW1u1b3azD4fdknBbEjacPfqp7jD2ox9iV4j6afufs7dfyfpfkmvCLym0H5pZs+XpMn/H61zozHizmUKJszMtLmveszdPxZ6PSG5+wfcfc3dD2jz38Q33b3ICc3dz0r6mZlVV/67XdL3Ay4ppNOSXmZmV0y+X25XoQ8uT/mqpLsnL98t6St1btTlqpC1tLxMQa5uk/RWSd81s+9M3vZBd/96wDUhDu+QdN9kAPqJpLcFXk8Q7v6gmR2W9JA2zy57WAVdhsDMviDpVZKeY2Ybkj4s6SOSvmRmb9fmnd9f1fpYXH4AAPLDM1QBIEPEHQAyRNwBIEPEHQAyRNwBIEPEHQAyRNwBIEP/Dz2jgNfaPCfLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAF41JREFUeJzt3WuI3Xedx/HP18REW9PUS7SaaZtapLMi67ZbdtWyItaCi2J94EIFpYhLWNj1DlJ9oM8WH4jog0XI1ttiUSQWFFd2DXVlWXYpG1tZq4mXaBJHExvZtQ1FE6vffTDnnzlzci7/++/2fkHpzGTOzG/+ybzPd37nf/5j7i4AQF6eEnoBAID+EXcAyBBxB4AMEXcAyBBxB4AMEXcAyNDKuJvZp83sUTN7ZOptzzKzI2b2o8n/nznsMgEATdSZ3D8r6bUzb7tH0gPu/iJJD0xeBwBEwuo8icnMDkj6mru/ZPL6DyS9yt3PmNnzJX3L3W8acqEAgPp2trzd89z9jCRNAv/cRe9oZgclHZQke+quP921b/u7+q50niFrFy30EgAU6MIvNn7l7vua3KZt3Gtz90OSDknS0/Zf69f/zXsv/dmF6y4O/el7t/v0rtBLAFCYH37ovaea3qZt3H9pZs+f2pZ5tMmNU4w6AKSk7amQX5V09+TluyV9pe4NCTsADK/OqZBfkPRfkm4ysw0ze7ukj0i6w8x+JOmOyesrpbS/DgApW7kt4+5vXvBHt/e8liRcuO4i++4AosczVAEgQ8QdADJE3FvgQWEAsSPuAJAh4t4S0zuAmBF3AMgQce+A6R1ArIg7AGSIuHfE9A4gRsS9BwQeQGyIe08IPICYEHcAyBBx7xHTO4BYEPeeEXgAMSDuAyDwAEIj7gMh8ABCGvwXZJeMX+zRvz2n0vhtXuevt9BLQOGI+8AIfDOpxHuVpl8HdwboG3EfAYGfL5eQ96HOseAOAE0Q95EQeGLe1bLjR/gxi7iPqMTAE/RxLDrORL9cxH1kpQSeqMeB6JeLuAeQc+CJehrm/T0R/LwQ90ByCjxBz8Ps3yOxTxtxDyj1wBP1vBH7tBH3wFINPGEvD7FPC3GPQHWpglQiT9ghbf93QOjjQ9wjksIUT9gxD6GPD3GPTMyBJ+yog9DHgbhHKMZtGsKONgh9OMQ9YrFM8YQdfSD04+J67pELfV14wo4h7Dnl/NsaGHFPQOjAA0Mh8sNhWyYRIfbhQ3zT7T1xYfTP+diNu0f/nNiOLZv+EffExLIP34cQIZ+nzjq4AxjPnlNO4HvQKe5m9h5Jfy3JJX1X0tvc/bd9LAyLjTHF9z21xxLytlatn/j3q/r3R+Tbax13M9sv6Z2SXuzuvzGzL0m6S9Jne1obVoh5ik895k0t+nqJfjdEvr2u2zI7JT3dzH4n6QpJv+i+JDQR2znxpUV9lXnHg+A3R+Sbax13d/+5mX1U0mlJv5H0DXf/xuz7mdlBSQclacezr2776bBCDFM8Ya+H4LfHfnx9rU+FNLNnSrpT0g2SXiDpSjN7y+z7ufshd7/V3W/dsefK9ivFSheuuxjktMm9Jy4Q9o6qY8ixXI3TJ+vpsi3zGkk/dfdzkmRm90t6haTP97EwtDfWVg0RGs7ssWWyvxxT/HJd4n5a0svM7AptbsvcLuloL6tCL4bcqiHs4yL287EXv1iXPfcHzeywpIckPSnpYUmH+loY+jHEFE/YwyP22zHFX67T2TLu/mFJH+5pLRhQ08ifv97Y10zIdOxLDT1T/HY8Q7UwTSI/L/DLpvZdxze6LU7SxfW1zh+jdKWHnil+E3EvVN3IL5vg+4h5049J/JspNfQEnrgXr+2e/BBh7/p5Cf9ypYW+9MATd0haHvnp6f2xG3dr74kLuri+Fizwi8xbD8Gfr5TQlxx44o5tZp8EVcV++hvk/PVP24z9jTdq74kL0UV+GsFfLffQlxp44o6lFj3j9cJ11Uu7tVfxTfHLEPzFqtDnFvkSA0/cccmBtXO13/eOa47ryNl1nfyLfZuTfAJT/DKz6y499jlO86UFnrhnpkmg67rjmuOXXn7TVQ9te/nwNbfoyNq6Tl63Gfm1xKb4RYj9lpym+ZICT9wTM0S8l1kU9kUeuzG9bZo6iH0+kS8l8MQ9YmOHfNp01KXFYa/efkTSSe2TtEvVPvw0Yp+PHCJfQuCJe2RCBl26POpNbT4Auxl4aSsEs/HLOfalhD71yOceeOIeWOiYS6uDXmc7Rtr8Wk5u7LsU+D2nfNs3/vSDdDGeJ9+X0qb61COfK+IeSApRb/qxjpxdvxR4aevc+OknQEmXT/O5Rr5SylSfYuRznt6J+8hSi3rdqX3W9LXk50W+lCl+VgmhTy3yuQa+9a/ZQ30H1s5d+i+0Pqf1SnUHUH3s6uucfQLU9DfQYzfu3vbNf3F9LdvYLbLr+Mal/3KU0q8MzPHy1kzuA4oh5pU2UW87tU+b/W1Qs1eZnDfFV3KN3jw5T/R7T1xIZorPCXEfSCxhbzupdw379N777EXJ5gVeuvxa8YQ+n9CnsFWT2/YM2zI9i2X7RQoT9unPuew4zPsmWvaNn0vkmspt6yb2rZqctmeY3HuUetSHNm+LRlKtKV4q5+yaRXKa6NmqGR6Tew9ymNYrbaf26dstm97nXWWyzRSfety6ymGij3WCz2V6J+4dxRT1UGFfpUvg2apZLeXIx7pNk0PgiXsHMYU9NXUDL7EXX1fK03yMgU8dcW8pt7D3PbXPrmve8eoz8ER+uxRDH9sUn/r0TtxbyC3sY6l73JYFnsg3l2Lk0R1xbyjHsA+1115njYt+jd+y841XnWVB5OdLKfKxBD7l6Z24N5Bj2PtU506i7vaM1C3wEvvxi6SyZRNL4FNF3BPTd9iHmtqbaht4pvhuYo98DPvwqU7vxL2mGKb2WCf2ZeatuemxXPWU8LpTPJFfLIXIoxniXkOuYY9laq8smt6lfgIvEflVYo58yMCnOL0T9xVyDfuYmkzvXQNP5PsRa+SZ4Osj7oUaampv8nGHCLzU7MqDRH65GCMfKvCpTe/EfQmm9nDGDLzEmTWrxBZ5JvjViHvEhgp7iL32RV9L2zvQoQJP5JeLKfIEfjnivkAMU3vplk3v0jCBl4h8HbFEfuzAp7Q10ynuZna1mR02s+NmdszMXt7XwkpX0nbMsjvSvgJP5IcRQ+SZ4OfrOrl/QtK/uPu6pJdKOtZ9SeHlPLXHdvrjmNr+cggCv1royI8Z+FSm99ZxN7OrJL1S0qckyd0vuvuv+1pYyVKf2hfdgbT9uvqY3isEflihI48tXSb3F0o6J+kzZvawmd1rZlfOvpOZHTSzo2Z29Pfnn+jw6cYRempPPextdT3uYwWeyNcTIvJsz2zXJe47Jd0i6ZPufrOkJyTdM/tO7n7I3W9191t37Lms/RhRylsyq6Z3aZzAS0zxTeQa+BS2ZrrEfUPShrs/OHn9sDZjDyy07CeTVdM7gU/T2FM8E/ym1nF397OSfmZmN03edLuk7/eyqkDYkilP18AT+frGjDyB7362zDsk3Wdm/yPpTyT9ffclAYv1Pb1L3QIvMcU3xQOu4+gUd3f/zmQ//Y/d/Y3u/n99Law0Q0/tqey39/XTU9PAd0Xgmxljih96eo99351nqE6E3pIpSdc7sjrTe1Ndp3eJwLeReuBjRtwxiC4/KYSa3gl8GJwbPwziHgEeSI0HgQ9nqMCXOr0TdySp7tbM2Hvv6Ca1CT7mfXfiXoBUHkyNBdN7WENs05Q4vRN38WBqCKu2our8ncQ+vRP4blKb4mND3IE5+pje0V2fgS9teifugfFg6jiY3tPFBN8OcQcW6Gt6J/DdEfjmiDuSNsQTmhCnPgJf0tYMcUe0eKAbKYj1dEjiDoyArZl+sD1TH3FHMdo8qMpZM/Eh8PUQdwDIEHEHkJwu03spD6oSdwDIEHEHgAwRdwDIEHEHkJwup5aWcgYUcQeADBF3AEnhCWH1EHcUo83TxPs6bY4n3vSDsNdH3BGtkxv7Qi8BWCnWX+VI3JG03ad3hV4CRtLH1F7Kg6kScQcWYksmHmzHNEfcUYRYL8uK1Qh7O8Q9sCNn10MvAXMwtcehz7CXtCUjEXcE0sedWt399lBTO2Hvhom9G+Ku/M/KOPz4LaGX0FjIv5NSrhoYs77DXtrULkk7Qy8AGBJTe1pSm9ZjPQ1SYnKPAvvuzQ11CmQfUzthb2eosJc4tUtM7hhIl62gvrZkmk7thD2M1Kb1VDC5T+S+7x6Trj+pxPrEJcLe3NBhL3Vql4h7NIbemknxQdUuxp7aCXszF9fXkg97zPvtEnFHZFb9BFVnaifscWMbZhyd425mO8zsYTP7Wh8LComtmfIQ9vGMMa1XSt6OqfTxgOq7JB2TdFUPH6toR86u645rjodexqCWbT+NPbV3CTtRr2/sSZ2wb+o0uZvZmqTXSbq3n+VgSCnvuxP2NOUa9tj326Xu2zIfl/R+SX9Y9A5mdtDMjprZ0d+ff6Ljpxte6K0Zznlvh7DHZcwtGMzXelvGzF4v6VF3/7aZvWrR+7n7IUmHJGn3DWtcmq8AbX5CWHan2uepj23DTtTrCRl0tmO26zK53ybpDWZ2UtIXJb3azD7fy6oCY3ofxlBfV92pnbAPJ/SkPmbYU9iSkTrE3d0/4O5r7n5A0l2Svunub+ltZYUbKoQx7rt3mdqHDPuu4xuEvYbQ2y9M7PNx+YEFTm7s04G1c6GXUbQ+wt426lgtdNSl8cOeytQu9fQkJnf/lru/vo+PhS05Te+Lvpa2W2CEPZzQWzAVJvbleIbqEqH33qV8999XWTa1DxF2tmBWiyXqEmGvg20Z9KrJTwWL7jzHDDtBXy2WoFdChT2lLRmJyX2lXKf3Mbdmmqy/S9j3nrhA2HsU06ReYWKvj8m9hhgeXM3t0gRN7zTrhL0uor5YbDGfFjLsqU3tEpN7Uvqe4MeY3uetuel2TF9hZ199sRin9GlM7M0R95pi2J6R8n2AdciwE/XFYo+6FD7sKU7tEnFvJMfAj31a5Lxj2CbsdfbXifpiqUQ9dNhTRtwbIvD1P06XNa4K+zJEfb4q6LFHXQo/rVdSndol4t5KjoEfQ92pvW3Yifp8qQS9EkvYU0fcW8ot8H1vz9RZV5Owr9qGIerbpTSlT4sp7ClP7RJx7yS3wA9p9lg1DfsiTOvbpRh0if31IRD3jmIKfNfID/XgatuwM63Xk+qUXokx6qlP7RJPYupFFa/QT3SSwj3ZafqOYdmdTJOwL0LUN6Ua80qMUZfyCLvE5N6rmKb4tvqe3qePSZ2wM60vl/qUXok17Dlhcu9ZDJcqkLYC32aKP/z4LXrTVQ91+rzS8js7pvX6Ug/5tNijnsvULhH3QcS2TSM1j3yXwEvL99nrhp2g54Wwj4u4Dyi2yI8deOnyrZh52zDTCHp+Yo96roj7CGKJfJetmmWqffrq41df77KwM6lvyjXoUlpRz21ql4j7qFKMfNvpvUnYS4t6zkGX0oq6lGfYJeIexPR+dMjQDzHJz+61E/X8Y15JLeq5I+6BxRD6VZGvM71PnyWz+/SuS1FfdAZM7mEvJeiVVMOe69QuEfeozE69Y8e+6yR/cmPfpbCXtqdeWswrqUZdyjvsEnGPWqjYz4v8oul99klPVdhzDnmFoKcr97BLxD0py54UNET4654+eeTs+qWpXfJsw15qzCs5RF0qI+wScc9G20sfrLpTmN5LnzfRHzm7rjP/uV/POeXae+K3WYW99JhXcom6VE7YJeJevDZ3Cv84c72YPadc+/75RJ/LCoKYb5dT1KWywi4Rd9Sw7HecpjqtE/L5cgt6pbSwS8Qdc6z6hdXVmTCpRJ2QL5dr0Cslhl0i7phYFHRp/hORYgw7EW8m96hL5YZdIu7FWxb1ZS6urwULPBFvr4SgV0oOu0Tci1U36ot+r6m0PbJ9h56A96ekoFdKD7tE3IvTZFKfF/bHbtw999mnxDguJQa9Qtg3EfdCNN1+WTaxIz4lx7xC1Lcj7plru6e+zKLpHeMi6FsI++Vax93MrpX0T5KukfQHSYfc/RN9LQzdDRH2CoEfHzGfj7DP12Vyf1LS+9z9ITPbI+nbZnbE3b/f09rQ0pBRn0bgh0XMlyPqy7WOu7ufkXRm8vJ5Mzsmab8k4h7QWGGvVAEi8t0R8/oI+2q97Lmb2QFJN0t6cM6fHZR0UJJ2PPvqPj4d5hg76rOIfDOEvB2iXl/nuJvZMyR9WdK73f3x2T9390OSDknS7hvWOAWjZ6GjPovIX46Qd0fUm+sUdzN7qjbDfp+739/PklBXbGGfNhu0EmJPxIdB2NvpcraMSfqUpGPu/rH+loQ6hg77+eut13Pdc4k9AR8PUe+my+R+m6S3SvqumX1n8rYPuvvXuy8Li8Q8rTcxL5Khg0+440DU+9HlbJn/kMTfwojGDnvf0/sqxLVsRL1fPEM1AblM68Asgj6cp4ReAJYLHXa++TCE89cb/7YGRtwjFjrsFb4J0ReiPh62ZSIUS9Snjb3/jnwQ8zCIe2RiDHuFwKMugh4ecY9IzGGvEHgsQtDjQtwjkULYKwQeEjGPHXEPLKWoTyPwZSLo6SDuaK36Rify+SLm6SLugaQ6sc9D5PNAyPNC3APIKezTiHw6CHn+iPvIcg37tOlwEPrwCHmZiPuISgj7LEI/HiKOacR9JCWGfdZsfIh9cwQcdRH3ERD2+Yj9doQbfSLuAyPs9S2LW8rhJ9oIgbgPiLD3p2sgm9w5EGPkgLgPhLDHhWCjNFzPfQCEHUBoxL1nhB1ADIh7jwg7gFgQdwDIEHHvCVM7gJgQdwDIEHHvAVM7gNgQ944IO4AYEfcOCDuAWBF3AMgQcW+JqR1AzIg7AGSIuLfA1A4gdsQdADJE3AEgQ8S9IbZkAKSAuANAhjrF3cxea2Y/MLMfm9k9K9//ojH5AsAIWsfdzHZI+gdJfynpxZLebGYvrnNbAg8Aw+oyuf+ZpB+7+0/c/aKkL0q6s+6Nd5/elVzkU1svgHJ1+QXZ+yX9bOr1DUl/PvtOZnZQ0sHJqxd++KH3PtLhc+bkOZJ+FXoRkeBYbOFYbOFYbLmp6Q26xH3er5P3y97gfkjSIUkys6PufmuHz5kNjsUWjsUWjsUWjsUWMzva9DZdtmU2JF079fqapF90+HgAgJ50ift/S3qRmd1gZrsk3SXpq/0sCwDQRettGXd/0sz+TtK/Stoh6dPu/r0VNzvU9vNliGOxhWOxhWOxhWOxpfGxMPfLtskBAInjGaoAkCHiDgAZGiXuTS9TkCszu9bM/s3MjpnZ98zsXaHXFJqZ7TCzh83sa6HXEpKZXW1mh83s+OTfx8tDrykUM3vP5PvjETP7gpk9LfSaxmJmnzazR83skam3PcvMjpjZjyb/f2adjzV43LtcpiBDT0p6n7v/kaSXSfrbgo9F5V2SjoVeRAQ+Ielf3H1d0ktV6DExs/2S3inpVnd/iTZP1rgr7KpG9VlJr5152z2SHnD3F0l6YPL6SmNM7p0uU5ATdz/j7g9NXj6vzW/g/WFXFY6ZrUl6naR7Q68lJDO7StIrJX1Kktz9orv/Ouyqgtop6elmtlPSFSro+TPu/u+S/nfmzXdK+tzk5c9JemOdjzVG3OddpqDYoFXM7ICkmyU9GHYlQX1c0vsl/SH0QgJ7oaRzkj4z2aK618yuDL2oENz955I+Kum0pDOSHnP3b4RdVXDPc/cz0uaAKOm5dW40RtxrXaagJGb2DElflvRud3889HpCMLPXS3rU3b8dei0R2CnpFkmfdPebJT2hmj9652ayn3ynpBskvUDSlWb2lrCrStMYcecyBVPM7KnaDPt97n5/6PUEdJukN5jZSW1u1b3azD4fdknBbEjacPfqp7jD2ox9iV4j6afufs7dfyfpfkmvCLym0H5pZs+XpMn/H61zozHizmUKJszMtLmveszdPxZ6PSG5+wfcfc3dD2jz38Q33b3ICc3dz0r6mZlVV/67XdL3Ay4ppNOSXmZmV0y+X25XoQ8uT/mqpLsnL98t6St1btTlqpC1tLxMQa5uk/RWSd81s+9M3vZBd/96wDUhDu+QdN9kAPqJpLcFXk8Q7v6gmR2W9JA2zy57WAVdhsDMviDpVZKeY2Ybkj4s6SOSvmRmb9fmnd9f1fpYXH4AAPLDM1QBIEPEHQAyRNwBIEPEHQAyRNwBIEPEHQAyRNwBIEP/Dz2jgNfaPCfLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "#Armamos una grilla de nxn llena de ceros\n",
    "\n",
    "dominio = np.zeros((n,n))\n",
    "curva_frontera1 = curva_frontera_tri1\n",
    "curva_frontera2 = curva_frontera_tri2\n",
    "\n",
    "ptsfrontera1 = [ curva_frontera1(t) for t in np.linspace(0.,1.,2 * n)]\n",
    "ptsfrontera2 = [ curva_frontera2(t) for t in np.linspace(0.,1.,2 * n)]\n",
    "valsfrontera1=[1. for p in ptsfrontera1]\n",
    "valsfrontera2=[-1. for p in ptsfrontera2]\n",
    "\n",
    "ptsfrontera = ptsfrontera1 + ptsfrontera2 \n",
    "valsfrontera = valsfrontera1 + valsfrontera2\n",
    "\n",
    "for p,v in zip(ptsfrontera,valsfrontera):\n",
    "    if 0<p[0]<n and 0<p[1]<n:\n",
    "        dominio[p[0],p[1]] = v\n",
    "    \n",
    "    \n",
    "X, Y = np.meshgrid(np.linspace(0,lado,n),np.linspace(0,lado,n))\n",
    "plt.contour(X, Y,dominio)\n",
    "\n",
    "\n",
    "def iterdom(dominio):\n",
    "    newdominio = dominio.copy()\n",
    "    for nx in range(1,n-1):\n",
    "        for ny in range(1,n-1):\n",
    "            newdominio[nx, ny] =  .25 * (dominio[nx + 1 ,ny] + \n",
    "                                      dominio[nx-1,ny] + \n",
    "                                      dominio[nx,ny+1] + \n",
    "                                      dominio[nx,ny-1])\n",
    "    for p in ptsfrontera:\n",
    "        if 0<p[0]<n and 0<p[1]<n:\n",
    "            newdominio[p[0], p[1]] =  dominio[p[0], p[1]]\n",
    "\n",
    "    err = (newdominio - dominio).flatten()\n",
    "    err = np.sqrt(err.dot(err))\n",
    "    return (newdominio,err)\n",
    "\n",
    "\n",
    "plt.contourf(X, Y, dominio)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "for i in range(100):\n",
    "    dominio,err = iterdom(dominio)\n",
    "\n",
    "plt.contourf(X, Y, dominio)\n",
    "plt.show()\n",
    "print(err)\n",
    "\n",
    "for i in range(5000):\n",
    "    dominio,err = iterdom(dominio)\n",
    "\n",
    "plt.contourf(X, Y, dominio)\n",
    "plt.show()\n",
    "print(err)\n",
    "\n",
    "for i in range(200):\n",
    "    dominio,err = iterdom(dominio)\n",
    "\n",
    "plt.contourf(X, Y, dominio)\n",
    "plt.savefig(\"equipot.svg\")\n",
    "plt.show()\n",
    "print(err)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Debido a que el promedio siempre se toma a primeros vecinos, partir de una solución con un potencial nulo salvo en los puntos que definen la condición de frontera, implica que la convergencia llevará una cantidad de pasos proporcional al tamaño de la grilla. Podemos mejorar esta convergencia si partimos de una mejor aproximación de la solución. \n",
    "\n",
    "Problema: Construir un algoritmo que resuelva el problema para un n=n0 pequeño (ej, n0=10), y a partir de esa solución, construya la grilla inicial para resolver el problema con un n1=2*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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}