tutorialIntegracionNumerica1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial de solución de ecuaciones diferenciales con parámetros constantes\n",
    "### Lourdes Martín-Aguilar & Marco Herrera Valdez\n",
    "#### Facultad de Ciencias, Universidad Nacional Autónoma de México\n",
    "\n",
    "Resolvamos un problema con valor inicial para una variable $x$ que depende de $t$ \n",
    "\\begin{equation}\n",
    "\\partial_t x = \\frac{a - x}{\\tau}, \\quad x(t_0) = 0 \n",
    "\\end{equation}\n",
    "utilizando las rutinas de integración numérica del módulo _scipy_ de _python_. \n",
    "\n",
    "La solución analítica es de la forma\n",
    "\\begin{equation}\n",
    "x(t) = a - \\left( a - x_0 \\right) \\exp \\left(-\\frac{t}{\\tau} \\right) \n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definiciones de funciones y simulación numérica\n",
    "\n",
    "El problema de valor inicial descrito arriba será resuelto mediante el uso de un iterador. Las herramientas necesarias para la integración numérica son parte del paquete _scipy_. La solución numérica y la analítica serán comparadas gráficamente. \n",
    "\n",
    "Primero importemos los módulos necesarios"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import scipy as sc\n",
    "import matplotlib.pylab as gr\n",
    "from scipy.integrate import odeint\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El módulo scipy tiene herramientas para hacer cálculos numéricos, el módulo pylab de matplotlib continene herramientas para graficar arreglos, listas, y otros objetos de python. El comando _matplotlib inline_ instruye al kernel de jupyter (o ipython notebook) para que muestre las figuras dentro del documento."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generemos un diccionario con los parámetros y arreglos necesarios para integrar la ecuación numéricamente."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['timeMax', 'ic', 'tau', 'a', 'sampTimes', 'timeStep'])\n",
      "[  0.00000000e+00   1.00000000e-01   2.00000000e-01 ...,   2.99700000e+02\n",
      "   2.99800000e+02   2.99900000e+02]\n"
     ]
    }
   ],
   "source": [
    "p=dict()\n",
    "p[\"a\"]=100.0; p[\"tau\"]=40.0; p[\"ic\"]=260.2\n",
    "p[\"timeMax\"]=300.0; p[\"timeStep\"]=0.1;\n",
    "p[\"sampTimes\"]= sc.arange(0,p[\"timeMax\"],p[\"timeStep\"])\n",
    "print(p.keys())\n",
    "print(p[\"sampTimes\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Definamos una función para que python la integre numéricamente"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def linearEq(u,t,p):\n",
    "    du = (p[\"a\"]-u)/p[\"tau\"]\n",
    "    return du"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "y otra función para la solución analítica"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f= lambda t,p: p[\"a\"] - (p[\"a\"]-p[\"ic\"])*sc.exp(-t/p[\"tau\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora iteremos la función usando la función _integrate.odeint_ de scipy. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "uSol= odeint(func=linearEq, y0=p[\"ic\"], t=p[\"sampTimes\"],args=(p,)).transpose()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El integrador da una solución en forma de un arreglo de $m \\times n$ donde $m$ es el número de puntos en el eje temporal considerados para la integración y $n$ es el número de variables de estado a integrar; en este caso sólo una. La transposición hace más fácil extraer las soluciones para las distintas variables de estado usando un sólo índice. \n",
    "\n",
    "Grafica la solución $y(t)$ y compara contra la solución numérica $u(t)$ y un acercamiento"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x103184f98>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApkAAAE4CAYAAAAQB+iEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VdW99/HPL0FGIYQxQhAOOBSvVRxKn9ZWghb1Wq0d\nnsttlVKubfW5PqIWayu23BCxDy3XYm17vW2t1lku9bZVbx0IQjpYAWWQIdQBEsQEEgIESMIU8nv+\nODsnJ+FkIJzknJN836/Xebn32nut89uI299Za6+1zd0REREREYmntEQHICIiIiJdj5JMEREREYk7\nJZkiIiIiEndKMkVEREQk7pRkioiIiEjcKckUERERkbhrNck0s2wzW2Zmm8xsg5ndFpRPNLNVZrY2\n+OfFUXVmm9l7ZrbZzK7oyAsQEelMMe6JM4PyXDP70MzWBJ+rgvJeZvaMma0P6tzdTLuZZrbEzN4x\ns1fNLKMzr0tEJN6stXUyzSwLyHL3dWZ2KvAW8AXgIWC+uy8xs38EvuPuk83sHOBp4GNANrAUONO1\nIKeIdAEx7omrgeuAfwYOuPvCJud/DbjS3a83sz5AITDJ3T9oct6PgN3uvsDMvgtkunvMhFREJBW0\n2pPp7jvdfV2wXQX8HRgB7AAGBqcNBEqC7c8Bi9y91t2LgfeAiXGOW0QkIWLcEzcDI4PDFqPKTqCf\nmaUDfYHDwP4Y510HPB5sPw58Pp5xi4h0thN6JtPMxgATgJXA3cCPzewDYAEwOzhtJLA9qloJDTdg\nEZEuo8k9EeBWM1tnZr82s4EA7v4q4aRyB1AM3O/ulTGaG+buZUGdncCwjo1eRKRjtTnJDIaFngNu\nD369PwLMdPfTgW8Bj3ZMiCIiySfGPfEhYKy7TyDce/nj4LxpQB8gCxgLfDtITlujR4xEJKX1aMtJ\nZtaD8M30SXd/Pij+uLtPAXD358zs10F5CTAqqno2DUPpEVdddZUfOnQosj9mzBhCodCJX0ECFBUV\npUys0VI1bkjd2FM1bkid2IuKiiguLo7s9+7dm1deeSXWsHXcxLonuvuuqFMeBl4Mtj8J/N7d64Bd\nZvY6cDHhXs1oZWY23N3Lguc+y2N9dyrfO9siVf7etZeuL3V1tWvrlHunu7f6AZ4AFjYpW0344XWA\ny4E3g+1zgLVATyAEvE8wwSj6M2nSJE9Vubm5iQ6hXVI1bvfUjT1V43ZP3diDe0ub7m3t/TRzT8yK\n2v4W8EywfRvwaLDdD9gEnBujzR8B3w22vwv8MNZ3p/K9sy1S9e9dW+n6UldXvjb3jrl3ttqTaWaX\nADcAG8xsLeEhnHuAm4CHzKwncCjYx90LzWwx4RmUR4Fb3DWzXES6hhbuideb2QSgjnAv5c1BlV8C\nj5jZBsITgx5x941BWw8D/+nuawgnmYvN7EZgGzC1865KRCT+Wk0y3f11IL2Zwx9vps58YP5JxCUi\nkpRauCe+0sz5h4FpzRz7ZtT2HuAz8YhRRCQZJOyNP2PGjEnUV5+0nJycRIfQLqkaN6Ru7KkaN6Ru\n7Kl8b2mLrn59qfr3rq10famrK18bdMy9JWFJZio/PJuqf9FSNW5I3dhTNW5I3dhT+d7SFl39+lL1\n711b6fpSV1e+NuiYe0ubZpeLSHz85Cc/Yd++fYkOo0vIyMjgjjvuSHQYItLBdN+Mr868dyrJFOlE\n+/btIzc3N9FhdAl5eXmJDkFEOoHum/HVmffOhA2Xi4iIiEjXpSRTREREROJOSaaIiIiIxJ2STBER\nERGJOyWZIiIiIhJ3SjJFpMMVFRXFLN+xYwc1NTWdHI2ISGpI9XunkkwRaeSee+7hpz/9adzaKyoq\nYuXKlTGPDR06lAULFjQq+/jHP87mzZvj9v0iIp1B987jKckUkYiKigqefPJJbr755ri1+Ytf/IIv\nf/nLMY/16NGDa665hieeeCJSdtdddzFnzpy4fb+ISEfTvTM2JZkiEvHYY49x9dVX06tXr7i0t379\nekaNGtXiORdffDFLly6N7F977bUsX76c8vLyuMQgItLRdO+MTUmmiES8/PLLTJo0qVHZ5ZdfTm1t\nbbvae/HFF5k8eXKr7Q0bNowtW7YA0KtXLy666CJeffXVdn2niEhn070zNr1WUiRJPPvss3Ft7ytf\n+coJ19mwYQNnn312ZL+kpAQID820x5tvvsk999zTanvnnXceq1evZty4cQCMHz+et99+m69+9avt\n+l4R6T5070zee6d6MkUkorKykv79+wOQn5/PrFmzyMrK4qmnnmq2TmlpKb///e8jN+a6ujpycnIA\nqKmpwcxabS8zM5MPP/wwst+/f38qKyvjeWkiIh1G987YlGSKSERmZiYHDhwAYMqUKaSnpzNr1iym\nTZvWbJ133nmHj33sY5SWlgLw1ltvMXr0aACOHTsWOa+l9vr06cORI0ci+wcOHGDgwIFxuy4RkY6k\ne2dsGi4XSRLtGaKJt/POO493332Xiy66CIB169ZFtpszefJkfvCDH3DDDTcA8Nprr3HFFVcAcMop\npzQ6t7n29u3bx6BBgyL7mzdvTprhHhFJbrp3Ju+9Uz2ZIhJx9dVXU1BQAEBhYSHjx48HYNGiRZFz\n3n//ferq6hrVW7lyJZ/61KcAWLp0KVOmTAEgKyuL6urqFtuD8MLC9c8UHT58mNWrV0faEBFJdrp3\nxqYkU0Qipk+fzssvv8zhw4cZNGgQGRkZLFq0KPKcEISXycjPz29U7/Of/zz/8z//w3/8x3+wZ88e\nhg0bBsCkSZNYtWoVQLPtQfhX+iWXXALACy+8wOTJk8nKyuq4CxURiSPdO2NTkikiEYMHD2b69On8\n8pe/JCsri0cffZQvf/nLjW5a69evjzzgDrBs2TLef/99vvOd71BRUcHtt98eOfaFL3yB1157DaDZ\n9g4fPkxGRga9e/cGYOHChdx7770dfakiInGje2dseiZTRBq57777Wjy+ZMkSrr766sj+4MGDOfvs\ns3n88ccJhUJMnz49cmzgwIEMHTqU3bt3M3jw4JjtPfvss9x0002R/TfeeOMkr0BEpPPp3nm8VpNM\nM8sGngCGA3XAw+7+0+DYTOAWoBb4o7vfHZTPBm4Mym939yUdE76IdLbPfvazjfbPP/98zj///GbP\nv+2223jkkUf4xje+cdyx7du3M2jQoEbry4mIdEXd8d7Zlp7MWmCWu68zs1OB1Wa2BMgCrgU+6u61\nZjYEwMzGA1OB8UA2sNTMznR3b9rwtGl5zJs3g1BodJwuR0SSjZnFvEkCjBo1qtVXp4mIdEdd4d7Z\n6jOZ7r7T3dcF21XAZmAk8K/AD929NjhWEVS5Dljk7rXuXgy8B0yM1fbTT3+bKVN+RlHRtpO+EBER\nERFJHic08cfMxgATgJXAWcClZrbCzJabWf0CTiOB7VHVSoKyGPqxZUsec+Y8dkJBi4iIiEhya3OS\nGQyVP0f4GcsqwkPtme7+v4DvAL9tXwj9KC2ta/00EZEkYGbZZrbMzDaZ2Ybg2XTMLNfMPjSzNcHn\nqqD8ejNbG5StNbNjZnZejHZj1hcRSVVtml1uZj0IJ5hPuvvzQfF24HcA7v5mcOMcTLjn8vSo6tlB\nWSNFRUXAXOAIFRV/o6Cg4Lj1n0REWlNQUBBZBBmguLi4o78y1nPq9YvfLXT3hdEnu/szwDMAZnYu\n8Ht3X99M28fVFxFJVW1dwuhRoNDdH4wq+wNwGfAnMzsL6Onuu83sBeBpM1tIeJj8DGBV0wZDoRDw\nbUKhOTz//OOa/CMi7ZKTk9PoB2peXl6Hfp+77wR2BttVZlb/nDqAtVL9K8CiFo63Vl9EJGW0Olxu\nZpcANwCXRQ35XAX8BhhrZhsI/0qfDuDuhcBioBB4Cbgl1sxygE9+8laeemqqEkwRSUlNnlMHuNXM\n1pnZr80sI0aVfwaebaHJ1uqLiKSMtswuf93d0919grtf4O4Xuvsr7n7U3b/q7h9194vd/U9Rdea7\n+xnuPr6lNTJvvfUKBg4cEK9rERHpNDGeU38IGOvuEwj3dC5scv5EoDr4IR5Li/VFRFJNwt/4c+DA\ngUSHICJyQmI9p+7uu6JOeRh4sUm1L9NCL2Yb6gPh59nnzp0b2W/6uICISFt0xvPsCU8y9+/fn+gQ\nRERO1HHPqZtZVvC8JsAXgY1Rx4zwSyo+1VyDLdWPFgqFyM3NPcnwRaS764zn2U9oncyOoCRTpHsI\nryhxvB07dlBTU9PJ0bRfC8+pLzCz9Wa2DpgEfCuq2qXAB8ELKqLbetjMLgx2W6ovIt1UKt87E55k\nHjlyhEOHDiU6DBEJ3HPPPfz0pz+Na5tFRUWsXLky5rGhQ4eyYMGCRmUf//jH2bx5c1xjiJcWnlOf\n7u7nBeWfd/eyqDp/cvdPxmjrm+6+Jthutr6IJD/dO4+X8CQT1JspkiwqKip48sknufnmm+Pa7i9+\n8Qu+/OUvxzzWo0cPrrnmGp544olI2V133cWcOXPiGoOISEfRvTM2JZkiSaCoaBvTpuUxeXIu06bl\nUVS0LSFtPPbYY1x99dX06tXrhOs2Z/369YwaNarFcy6++GKWLl0a2b/22mtZvnw55eXlcYtDRKSj\n6N4Zm5JMkQQrKtrGlCk/4+mnv01BQR5PP/1tpkz52QklifFoA+Dll19m0qRJjcouv/xyamtrT6id\naC+++CKTJ09utc1hw4axZcsWAHr16sVFF13Eq6++2u7vFRHpLLp3xqYkUyTB5sx5jC1b8oB+QUk/\ntmzJY86cxzq1DYANGzZw9tlnR/ZLSsJvhO3Ro/0LUbz55pucc845rbZ53nnnsXr16sj++PHjefvt\nt9v9vSIinUX3ztiSIsnUWpnSnZWU1NGQHNbrR2lpXae2AVBZWUn//v0ByM/PZ9asWWRlZfHUU081\nW2f79u387ne/4ytf+QoAtbW1fOYzn4kcP3jwIOEVfFpuMzMzkw8//DCy379/fyorK08ofhHpnszi\n82kv3TtjS2iSWf+HV1VVxbFjxxIZikjCjByZBlQ3Ka1mxIi2/+cZjzYgfLOq/9E3ZcoU0tPTmTVr\nFtOmTWu2zjvvvMPEiRMpLS0FYOXKlYwZMyZyPPq/7Zba7NOnD0eOHInsHzhwgIEDB55Q/CLSPbnH\n59NeunfGltAks1+/hp4X9WZKdzVv3gzGjculIUmsZty4XObNm9GpbUB42OXdd9+N7K9bt46LLrqo\nxTqf+cxnePzxx7nhhhsAeO211xr9Gm86tNNcm/v27WPQoEGR/c2bN3P++eefUPwiIomge2dsCU0y\nBwxoeG+5nsuU7ioUGk1+/kxuuOF+Jk/O5YYb7ic/fyah0OhObQPg6quvjrxmrLCwkPHjxwOwaNGi\nyDnvv/8+dXWNh+HfeOMNPvWp8Mts8vPzueyyy1iyZAkAWVlZVFdXt9gmhBcWHjduHACHDx9m9erV\nTJky5YTiFxFJBN07Y0voayUHDBgQ6SZWkindWSg0mqeeOrlXBcajjenTp3PBBRdw+PBhBg0aREZG\nBosWLWr06rFrr72Wn/zkJ1x55ZWRsi984Qu8+OKLbNq0iXHjxvHSSy9x1VVXATBp0iRWrVrF5MmT\nm20Twr/Sv/GNbwDwwgsvMHnyZLKysk7qekREOoPunbElPMmspyRTJPEGDx7M9OnT+eUvf8ltt93G\no48+etw569ev580332xU9vWvfz2y/U//9E+Njn3hC1/g/vvvj9z4YrV5+PBhMjIy6N27NwALFy7k\nkUceiccliYh0ON07Y0uaJHPfvn0JjERE6t13330tHl+yZAlXX311m9sbOHAgQ4cOZffu3QwePDjm\nOc8++yw33XRTZP+NN95oc/siIslA987jJfSZzIyMjMj2/v378ZOZ2iUineKzn/1sZGWItrrtttv4\n/e9/H/PY9u3bGTRoUKM15kREuprueO9MaE9mz5496dOnDwcPHqSuro6qqqrIOlMi0nWYWeSZoaZG\njRrV6qvTRES6o1S/dyZ8Mfbo3kwNmYuIiIh0DUoyRURERCTulGSKiIiISNwpyRQRERGRuEuqJHP/\n/v3HrYYvIiIiIqmn1dnlZpYNPAEMB+qAh939p1HH7wT+HRji7nuCstnAjUAtcLu7L2mu/VNOOYW+\nfftSU1MTmWEevX6mSFeSkZFBXl5eosPoEqJ/oIpI16X7Znx15r2zLUsY1QKz3H2dmZ0KrDazJe7+\n9yABnQJsqz/ZzMYDU4HxQDaw1MzO9BYWwczIyKCmpgYID5kryZSu6o477kh0CCIiKUX3zdTV6nC5\nu+9093XBdhWwGRgZHH4AuKtJleuARe5e6+7FwHvAxJa+Q89lioiIiHQtJ/RMppmNASYAK83sc8B2\nd9/Q5LSRwPao/RIaktKYlGSKiIiIdC1tfuNPMFT+HHA7cAy4h/BQ+UlTkikiIiLStbQpyTSzHoQT\nzCfd/XkzOxcYA7xt4RdxZgNrzGwi4Z7L06OqZwdljRQVFTF37lwA6urqcHfOOeecyAzztLSET3wX\nkRRQUFBAQUFBZL+4uDhhsYiISIO29mQ+ChS6+4MA7r4RyKo/aGZFwIXuvtfMXgCeNrOFhIfJzwBW\nNW0wFAqRm5sb2X/hhReorq7G3Tlw4IBmjopIm+Tk5JCTkxPZ1yxUEZHk0Gp3oZldAtwAXGZma81s\njZld1eQ0BwzA3QuBxUAh8BJwS0szy+tpyFxERESk62jL7PLX3T3d3Se4+wXufqG7v9LknLH1a2QG\n+/Pd/Qx3H9/SGpnRlGSKSCows2wzW2Zmm8xsg5nNDMpzzezD4Id45Me4mV0f9QN9rZkdM7PzYrSb\naWZLzOwdM3vVzDScIyIpLWkefIxOMisrKxMYiYhIi+rXDv4H4BPArWb2keDYwuCHeOTHuLs/U/8D\nHfgqsNXd18do925gqbufDSwDZnf8pYiIdJykSTIHDhwY2VaSKSLJqpW1g62V6l8BFjVz7Drg8WD7\nceDzJxmqiEhCJU2SOWDAAMIT1aGqqora2toERyQi0rLotYODolvNbJ2Z/bqZ4e5/Bp5tprlh7l4G\n4UQWGBbncEVEOlXSJJnp6emNXiep3kwRSWbRawcHPZoPAWPdfQKwE1jY5PyJQHUwObItWp0wKSKS\nzNq8GHtnGDhwYGTST2VlJUOGDElwRCIix2u6djCAu++KOuVh4MUm1b5M872YAGVmNtzdy8wsCyiP\ndVL0GsNw/BJOIiJt0RlrDCdVkpmZmcm2bdsA9WSKSFJrtHYwgJllBcPcAF8ENkYdM2Aq8KkW2nwB\nmAH8CPga8Hysk5quMSwi0h6dscZwUiWZ0ZN/9u7dm8BIRERii1o7eIOZrSU8rH0PcL2ZTQDqgGLg\n5qhqlwIfuHtxk7YeBv7T3dcQTi4Xm9mNwDbCSamISMpK2iSzsrISd49MBhIRSQbu/jqQHuPQKzHK\n6uv8CfhkjPJvRm3vAT4TjxhFRJJB0kz8AejTpw+9evUCoLa2lurq6gRHJCIiIiLtkVRJJmi9TBER\nEZGuIOmSzMzMzMi2nssUERERSU1Jl2SqJ1NEREQk9SVdkhndk6kkU0RERCQ1JV2S2fT1kkePHk1w\nRCIiIiJyopIuyUxLSyMjo+GVv+rNFBEREUk9SZdkgp7LFBEREUl1SbUYe73MzExWrXqTxYvXcujQ\nKsaPH8i8eTMIhUYnOjQRERERaYOkTDL37z/A/PmbKCv7OdCPNWuqWbEil/z8mUo0RURERFJAUg6X\nP/DAi5SVLQT6BSX92LIljzlzHktgVCIiIiLSVkmZZO7caTQkmPX6UVpal4hwREREROQEJWWSOXJk\nGtD0veXVjBiRlOGKiIiISBNJmbXNmzeDUaPupiHRrGbcuFzmzZuRsJhEREREpO1aTTLNLNvMlpnZ\nJjPbYGYzg/IFZrbZzNaZ2X+b2YCoOrPN7L3g+BUnGlQoNJrf/nY6n/zkrZxzzs1cdtmdmvQjIiIi\nkkLaMru8Fpjl7uvM7FRgtZnlA0uAu929zsx+CMwGZpvZOcBUYDyQDSw1szPd3U8ksPPP/yi33hrO\nT9PT0xkz5vQTqS4iIiIiCdRqT6a773T3dcF2FbAZGOnuS929fibOCsIJJcDngEXuXuvuxcB7wMQT\nDax379707dsXgGPHjrFv374TbUJEREREEuSEnsk0szHABGBlk0M3Ai8F2yOB7VHHSoKyEzZo0KDI\n9t69e9vThIiIiIgkQJuTzGCo/Dng9qBHs778e8BRd3823sFlZmZGtvfs2RPv5kVERESkg7TpjT9m\n1oNwgvmkuz8fVT4DuBq4LOr0EmBU1H52UNZIUVERc+fOjezn5OSQk5PT6Bz1ZIpIawoKCigoKIjs\nFxcXJyyWztKrF6SlhT/p6Q3bre2fyLknU1ffE3vbLPwR6S7a+lrJR4FCd3+wvsDMrgLuAi5198NR\n574APG1mDxAeJj8DWNW0wVAoRG5ubotfGp1k7tmzB3fH9F+oiERp+gM1Ly8vccF0kv37oa4Ojh0L\n/7P+09L+iZx7MnVP9nuOHu1a1xN9zL1xMpqMiXB3+54TrZuWluj/+lNLq0mmmV0C3ABsMLO1gAPf\nA34K9ATyg8Rvhbvf4u6FZrYYKASOArec6Mzyer1796ZPnz4cPHiQY8eOsX//fjIyMtrTlIhIl9Gr\nV6IjkPZwT/5E+GTPra1NnetpT13ougl3z57x/zvfapLp7q8D6TEOndlCnfnA/JOIK2LQoEGUlIRH\n2/fu3askU0REUpJZ+H/q6bH+jyopof6HQjIl9vH6nvokOp7aOlyeMJmZmZEkc8+ePYwZMyaxAYmI\niEi31JV/KHTEk0ZJ/3RB9HOZu3fvTmAkIiIiItJWSZ9kDh48OLK9d+9e6urqEhiNiIiIiLRF0ieZ\nvXv3pl+/foDe/CMiiWdm2Wa2zMw2mdkGM5sZlOea2Ydmtib4XBVV5zwz+5uZbTSzt83suEfsW6ov\nIpKKkv6ZTAj3ZlZXVwNQUVHRaJF2EZFOVgvMcvd1wUsqVptZfnBsobsvjD7ZzNKBJ4Eb3H2jmWUS\nXnkjluPqi4ikqqTvyYTGQ+Z6LlNEEsndd7r7umC7CthMw6tzYy3kewXwtrtvDOrsbWFZNy0ELCJd\nhpJMEZF2MrMxwARgZVB0q5mtM7Nfm1n9emtnBee+YmZvmdldLTQZq76ISEpKieHyQYMGYWa4O/v3\n7+fo0aOccsopiQ5LRLqxYKj8OeB2d68ys4eAe93dzew+YCHwdcL32UuAi4FDwGtm9pa7L2/SZHP1\nG2nLK3lFRFrTGa/kTYkkMz09nYEDB0beX757926ysrISHJWIdFdm1oNwgvmkuz8P4O67ok55GHgx\n2P4Q+LO77w3qvgRcCDRKMluo30hbXskrItKapj9QO+KVvCkxXA4aMheRpPIoUOjuD9YXmFn0L98v\nAhuD7VeBj5pZ7yA5nUT4tbuNtFBfRCQlpURPJoSTzPfffx9QkikiiWNmlwA3ABvMbC3gwD3A9WY2\nAagDioGbAdy90swWAm8Fx15y95eDth4G/tPd1wALYtUXEUlVKZNkDhkyJLKtJFNEEsXdXwdivVTu\nlRbqPAM8E6P8m1Hb0+MSoIhIkkiZ4fL+/ftHJvscOnQosm6miIiIiCSflEkyzUzPZYqIiIikiJRJ\nMqHx5J+KiooERiIiIiIiLVGSKSIiIiJxl1JJZvTkn71793Ls2LEERiMiIiIizUmpJLNXr14MGDAA\ngLq6Oj2XKSIiIpKkUirJBBg6dGhke9euXS2cKSIiIiKJktJJpp7LFBEREUlOKZdkRj+XWVFRgbsn\nMBoRERERiSXlksz+/fvTu3dvAI4cOcL+/fsTHJGIiIiINNVqkmlm2Wa2zMw2mdkGM7stKM80syVm\n9o6ZvWpmGVF1ZpvZe2a22cyuiHfQei5TREREJLm15d3ltcAsd19nZqcCq81sCfAvwFJ3X2Bm3wVm\nA3eb2TnAVGA8kA0sNbMzPY7j2kOGDGH79u1AOMk844wz4tW0iEjSKysrw8xifoBmy9PS0k74uIhI\ne7WaZLr7TmBnsF1lZpsJJ4/XAZOC0x4HCoC7gc8Bi9y9Fig2s/eAicDKeAWtnkwR6c6WLVvWqd/X\nXELbUlKb6OPQvqQ60XGfyHGRZNeWnswIMxsDTABWAMPdvQzCiaiZDQtOGwm8EVWtJCiLm8zMTNLT\n0zl27BjV1dUcPHiQPn36xPMrREQk4O6aZJmEultSnSrHpUGbk0wLD5U/B9we9Gg2veN02h0oLS2N\nwYMHU15eDoR7M08//fTO+noRkYQaNmxYJOmrTwBjfeJxXJJX039Pegtecki2pLetxztCm5JMM+tB\nOMF80t2fD4rLzGy4u5eZWRZQHpSXAKOiqmcHZY0UFRUxd+7cyH5OTg45OTltDnzo0KFKMkWEgoIC\nCgoKIvvFxcUJi6WzXH755Z32XU2Tzrq6uuPKT+R4RyfF3eW4JK/of1fdXVt7Mh8FCt39waiyF4AZ\nwI+ArwHPR5U/bWYPEB4mPwNY1bTBUChEbm5uO8OuTzJ3sXjxWqqqXuejH/0f5s2bQSg0ut1tikjq\nafoDNS8vL3HBdEFNeznS09MTGI3US7aktzOPw8n/2OmoH0PSWKtJppldAtwAbDCztYSHxe8hnFwu\nNrMbgW2EZ5Tj7oVmthgoBI4Ct3gH/MlXVdUwf/4mysp+DvRj/fpqVqzIJT9/phJNERHp0qKHPCW5\npGpS/e6778b9z6Its8tfB5r76fqZZurMB+afRFytyst7irKyhUC/oKQfW7bkMWfO/Tz1VPt7SEVE\nRETaSxOAGqTcG3/qlZTU0ZBg1utHaWldIsIRERERkSgpm2SOHJkGVDcprWbEiJS9JBEREZEuI2Uz\nsnnzZjB27L/RkGhWM3bsHObNm5GwmEREREQk7IQWY08modBoli69jW9840527nQGDjzCggU3adKP\niIiISBJI2SQTwonmAw/cwqZNmwDo1euUBEckIiIiIpDCw+X1hg0bFtmuX5xdRERERBIr5ZPMoUOH\nRt7Hum/fPg4dOpTgiEREREQk5ZPM9PR0Bg8eHNlXb6aIdCQzyzazZWa2ycw2mNnMoDzXzD40szXB\n56qoOufQD/RFAAAgAElEQVSZ2d/MbKOZvW1mPWO0m2lmS8zsHTN71cwyOvO6RETiLeWTTGg8ZF5W\nVpbASESkG6gFZrn7PwCfAG41s48Exxa6+4XB5xUAM0sHngRucvdzgRzCb0Nr6m5gqbufDSwDZnfw\ndYiIdKgukWRmZWVFtnfu3JnASESkq3P3ne6+LtiuAjYDI4PDsV7zcQXwtrtvDOrsbeZVu9cBjwfb\njwOfj2vgIiKdrEskmUOGDKFHj/BE+aqqKqqqqhIckYh0B2Y2BpgArAyKbjWzdWb266jh7rOCc18x\ns7fM7K5mmhvm7mUQTmSBYc2cJyKSElJ6CaN6aWlpDBs2jNLSUiDcm3nGGWckOCoR6crM7FTgOeB2\nd68ys4eAe93dzew+YCHwdcL32UuAi4FDwGtm9pa7L2/lK2L1dlJUVMTcuXMj+zk5OeTk5Jzs5YhI\nN1NQUEBBQUFkv7i4OO7f0SWSTAgPmSvJFJHOYGY9CCeYT7r78wDuvivqlIeBF4PtD4E/u/veoO5L\nwIVA0ySzzMyGu3uZmWUBMWcxhkIhcnNz43cxItItNf2BmpeXF/fv6BLD5QCnnXZaZHvnzp3EfuRJ\nRCQuHgUK3f3B+oIgMaz3RWBjsP0q8FEz6x0kp5OAwhhtvgDMCLa/Bjwf76BFRDpTl0kyBwwYQN++\nfQE4evQou3fvTnBEItIVmdklwA3AZWa2Nmq5ogVmtt7M1hFOJL8F4O6VhIfO3wLWAKvd/eWgrYfN\n7MKg6R8BU8zsHeBy4IedemEiInHWZYbLITxkvnXrViDcmzlkyJAERyQiXY27vw6kxzj0Sgt1ngGe\niVH+zajtPcBn4hGjiEgy6DI9maCljERERESSRZdNMisqKjh6NNZ6xyIiIiLS0bpUktmrVy8yMzMB\ncHe9/UdEREQkQbpUkgmNZ5nv2LEjgZGIiIiIdF9dOsmsXzdTRERERDpXl0syhw4dSs+ePQGoqamh\nsrIywRGJiIiIdD+tJplm9oiZlZnZ+qiyj5nZqmCNuFVmdnHUsdlm9p6ZbTazKzoq8BbibdSbWVJS\n0tkhiIiIiHR7benJ/A1wZZOyBcD33f0CIBf4dwAzOweYCowH/hF4yMwsfuG2zYgRIyLbGjIXERER\n6XytJpnu/ldgb5PiHUBGsD0QqO8u/BywyN1r3b0YeA+YGJ9Q2+60006jPretqKjg8OHDnR2CiIiI\nSLfW3mcy7wYWmtkHhHs1ZwflI4HtUeeVBGWdqlevXgwePDiyr1nmIiIiIp2rvUnmI8BMdz+d8Pt5\nH41fSPExcmRDbqvnMkVEREQ6V3vfXf5xd58C4O7Pmdmvg/ISYFTUedk0DKU3UlRUxNy5cyP7OTk5\n5OTktDOc440YMYK3334bCL9isq6ujrS0LjeZXqTbKygooKCgILJfXFycsFhERKRBW5NMCz713jOz\nSe7+JzO7nPCzlwAvAE+b2QOEh8nPAFbFajAUCpGbm9vOsFs3cOBA+vbtS01NDUeOHKGiooJhw4Z1\n2PeJSGI0/YGal5eXuGBERCSi1STTzJ4BcoDBwTOYucBNhGeO9wQOBfu4e6GZLQYKgaPALe7uHRR7\nq0aMGMH7778PhIfMlWSKiIiIdI5Wk0x3v76ZQx9v5vz5wPyTCSpesrOzI0nmhx9+yAUXXJDgiERE\nRES6h/Y+k5kShg8fzp49lTzzzCoqK3vxxBMrWLDgJkKh0YkOTURERKRL69JJ5rZt25k/fxMlJT8H\n+lFYWM3atbnk589UoikiIiLSgbr0dOs5cx6jpORHQL+gpB9btuQxZ85jCYxKREREpOvr0klmSUkd\nDQlmvX6UltYlIhwRERGRbqNLJ5kjR6YB1U1KqxkxoktftoiIiEjCdelsa968GYwbl0tDolnNaafd\nxbx5MxIWk4iIiEh30KUn/oRCo8nPn8k99/yI9evLGTjwMFOnXsDw4UMTHZqIiIhIl9alk0wIJ5rP\nPnsvy5cvZ+fOnUB4zcyzzjorwZGJiIiIdF1derg82qhRDa9U37ZtWwIjEREREen6ulWSaRZ+/XpF\nRQXV1U0nBImIiIhIvHSbJLNXr15kZWVF9j/44IMERiMiIiLStXWbJBNg9OiGt/xoyFxE2sPMss1s\nmZltMrMNZjYzKM81sw/NbE3wuSooH21mNVHlDzXTbsz6IiKpqstP/ImWnZ1NWloadXV17N27l/37\n9zNgwIBEhyUiqaUWmOXu68zsVGC1meUHxxa6+8IYdd539wvb0HZz9UVEUk636sk85ZRTGDlyZGRf\nvZkicqLcfae7rwu2q4DNQP2NxZqp1lx5e88TEUl63SrJhMZD5nouU0ROhpmNASYAK4OiW81snZn9\n2swGRp06JhgCX25mn2qhyej6GR0UtohIp+hWw+UAI0aMoEePHtTW1rJ//3727t1LZmZmosMSkRQT\nDJU/B9zu7lXBs5b3urub2X3Aj4GvAzuA0919r5ldCPzBzM4JekGjNa2/MKjfSFFREXPnzo3s5+Tk\nkJOT0wFXKCJdWUFBAQUFBZH94uLiuH9Ht0sy09PTyc7OjvxhFhUVKckUkRNiZj0IJ5hPuvvzAO6+\nK+qUh4EXg/IjwJFge42ZbQHOAtZEt9lc/aZCoRC5ublxuhIR6a6a/kDNy8uL+3d0u+FygDFjxkS2\nt23bRl1dXeKCEZFU9ChQ6O4P1heYWVbU8S8CG4PyIWaWFmyPBc4AtjZtsLn6IiKpqtv1ZAJkZWXR\np08fDh48yKFDh9ixY0ejCUEiIs0xs0uAG4ANZrYWcOAe4HozmwDUAcXAzUGVS4F7zexIcOxmd68M\n2noY+E93XwMsaKa+iEhK6pZJppkRCoUoLCwEwkPmSjJFpC3c/XUgPcahV5o5/3fA75o59s2o7elx\nCVBEJEl0y+FyCD/XVK+kpITDhw8nMBoRERGRrqXbJpkDBgxgyJAhANTV1WnNTBEREZE4ajXJNLNH\nzKzMzNY3KZ9pZpuD16r9MKp8tpm9Fxy7oiOCjpfo3syioqIERiIiIiLStbSlJ/M3wJXRBWaWA1wL\nfNTdPwrcH5SPB6YC44F/BB4ys6R9g8Xpp59OWlr4j2DPnj1UVlYmOCIRERGRrqHVJNPd/wrsbVL8\nr8AP3b02OKciKL8OWOTute5eDLwHTIxfuPHVs2dPRo0aFdnfsmVLAqMRERER6Tra+0zmWcClZrYi\neE3aRUH5SGB71HklNLzTNymNHTsWgPLyXdxxx0Pk5Pwb06blUVSkZzRFRERE2qu9Sxj1ADLd/X+Z\n2ceA3wJj4xdW5xk+fDjV1QeZP38TZWU/BfoB1axYkUt+/kxCodGtNSEiIiIiTbQ3ydxOsO6bu79p\nZsfMbDDhnsvTo87LDsqOkyzv3zUznn9+I2VlCwknmAD92LIljzlz7uepp/T6NpFk1hnv3xURkRPX\n1iTTgk+9PwCXAX8ys7OAnu6+28xeAJ42s4WEh8nPAFbFajCZ3r+7f38/GhLMev0oLdXrJkWSXWe8\nf1dERE5cq0mmmT0D5ACDzewDIJfwe3t/Y2YbgMPAdAB3LzSzxUAhcBS4xd29g2KPm1GjegDVNE40\nqxkxotsuIyoiIiJyUlpNMt39+mYOfbWZ8+cD808mqM42b94MXn/9+xQX30f9M5ljx/4b8+bdluDI\nRERERFJTt3x3eVOh0GiWLbuDG2/8FuXlxsCBR7j33uma9CMiIiLSTkoyA6HQaH71q7t46623ADh4\nsBp3J4nXkhcRERFJWnroMEooFKJnz54AHDhwgNLS0gRHJCIiIpKalGRG6dGjB+PGjYvsv/POOwmM\nRkRERCR1Kcls4qyzzooMkZeVlel95iIiIiLtoCSzib59+zZ6n7l6M0VEREROnJLMGM4+++zIdnFx\nMYcOHUpgNCIiIiKpR0lmDEOGDGHIkCEA1NXV8e677yY4IhEREZHUoiSzGR/5yEci2++++y5Hjx5N\nYDQiIiIiqUVJZjOys7MZMGAAAEePHuW9995LcEQiIiIiqUNJZjPMjPHjx0f2//73v3Ps2LEERiQi\nIiKSOpRktmDMmDH07dsXgMOHD7Nly5YERyQiIiKSGpRktiAtLa1Rb+bmzZupq6tLYEQiIiIiqUFJ\nZivGjRtH7969AaipqaGoqCjBEYmIiIgkPyWZrUhPT2+0bubGjRvVmynSjZlZtpktM7NNZrbBzGYG\n5blm9qGZrQk+VwXlo82sJqr8oWbazTSzJWb2jpm9amYZnXldIiLxpiSzDc466yx69epFefkuFiz4\nA5/4xHeZNi2PoqJtiQ5NRDpfLTDL3f8B+ARwq5nVr3m20N0vDD6vRNV5P6r8lmbavRtY6u5nA8uA\n2R12BSIinaBHogNIBT169GDAgIHMn/8qZWU/B/qxalU1K1bkkp8/k1BodKJDFJFO4u47gZ3BdpWZ\nbQZGBoetmWrNlUe7DpgUbD8OFBBOPEVEUpJ6Mtvo0Uf/RFnZQqBfUNKPLVvymDPnsQRGJSKJZGZj\ngAnAyqDoVjNbZ2a/NrOBUaeOCYbKl5vZp5ppbpi7l0EkkR3WUXGLiHQGJZlttGMHNCSY9fpRWqrn\nM0W6IzM7FXgOuN3dq4CHgLHuPoFwT+ePg1N3AKe7+4XAncAzQd3WeAeELSLSaTRc3kYjR6YB1TRO\nNKsZMUJ5ukh3Y2Y9CCeYT7r78wDuvivqlIeBF4PyI8CRYHuNmW0BzgLWNGm2zMyGu3uZmWUB5bG+\nu6ioiLlz50b2c3JyyMnJicdliUg3UlBQQEFBQWS/uLg47t+hJLON5s2bwYoVuWzZkkc40axm+PA7\nueeeOxIcmYgkwKNAobs/WF9gZlnBMDfAF4GNQfkQYI+715nZWOAMYGuMNl8AZgA/Ar4GPB/ri0Oh\nELm5ufG6DhHpppr+QM3Ly4v7dyjJbKNQaDT5+TP5/vf/nY0bKzj11GqmTr2Affv2Jjo0EelEZnYJ\ncAOwwczWEh7Wvge43swmAHVAMXBzUOVS4F4zOxIcu9ndK4O2Hgb+093XEE4uF5vZjcA2YGrnXZWI\nSPy1mmSa2SPANUCZu5/X5NidwL8DQ9x9T1A2G7iR8DIft7v7krhHnSCh0Giefnou5eXlvPbaa0C4\ne/nss89m0KBBCY5ORDqDu78OpMc49EqMMtz9d8Dvmjn2zajtPcBn4hGjiEgyaMsDhb8BrmxaaGbZ\nwBTCv7jry8YT/vU9HvhH4CEza8vSHSll2LBhZGdnR/bXrGn6aJWIiIhI99ZqkunufwVijQk/ANzV\npOw6YJG717p7MfAeMPFkg0xGEyZMIC0t/Me3a9cutm/fnuCIRERERJJHu6ZGm9nngO3uvqHJoZFA\ndLZVQsMixV1K//79OfPMMyP7a9asoba2NoERiYiIiCSPE04yzawP4Yfcu/30xnPPPZdevXoBUFNT\nw8aNGxMckYiIiEhyaM/s8nHAGODt4HnLbGCNmU0k3HN5etS52UHZcbrCWm89e/bkggsuYMWKFQD8\n/e9/Z8yYMQwcOLCVmiISL52x1puIiJy4tiaZFnxw941AVuSAWRFwobvvNbMXgKfNbCHhYfIzgFWx\nGuwqa72FQiG2bt1KeXk57s5bb73F5ZdfThec7ySSlDpjrTcRETlxrQ6Xm9kzwN+As8zsAzP7lyan\nOA0JaCGwGCgEXgJucfcu/2q0iy++OJJU7tq1i6KiogRHJCIiIpJYrfZkuvv1rRwf22R/PjD/JONK\nKRkZGYwfP57CwkLKy3cxY8YPMDuNUaNOYd68GYRCoxMdooiIiEin0ht/4uTcc89l1aq3mD9/E2Vl\nP6X+1ZMrVuSSnz9TiaaIiIh0K+1awkiOl56eziuvvE9Z2ULCCSZAP7ZsyWPOnMcSGJmIiIhI51OS\nGUd79vSkIcGs14/S0rpEhCMiIiKSMEoy42jkyDSguklpNSNG6I9ZREREuhdlP3E0b94Mxo3LpSHR\nrGb48Fn8y79cmsCoRERERDqfJv7EUSg0mvz8mcyZcz/vvLOPnj33MnXqBVRUlFNZWalF2kVERKTb\nUJIZZ6HQaJ56Kpe6ujqWLFnC3r17qaur4/XXX+fKK6+kRw/9kYuIiEjXp+HyDpKWlsYnP/lJ0tPT\nAdi/fz9r1qxJcFQiIiIinUNJZgcaMGAAF198cWR/y5YtehuQiIiIdAtKMjvY2LFjGT26YSH2VatW\nsWfPngRGJCIiItLxlGR2gokTJ5KRkQFAXV0df/nLXzh8+HCCoxIRERHpOEoyO0GPHj349Kc/zSmn\nnAJATU0Nr7/+Ou6e4MhEREREOoaSzE7Sv39/PvGJTwBQXr6LOXOe5GMfm8W0aXkUFW1LcHQiIiIi\n8aX1dDrRyJEjycwczLe+VUBZ2c+BfqxeXc2KFbnk588kFBrdahsiIiIiqUA9mZ3sySf/RlnZQhre\ncd6PLVvymDPnsQRGJSIiIhJfSjI7WWmp05Bg1utHcfGhRIQjIiIi0iGUZHaykSPTaHi3eb1q0tLC\nr54UERER6QqUZHayefNmMG5cLg2JZjXDh8/iS186n+XLl7N///4ERiciIiISH0oyO1koNJr8/Jnc\ncMP9TJ6cy//+3z9gzpwJDBs2lEOHDrF8+XKqq5v2dIpIsjCzbDNbZmabzGyDmc0MynPN7EMzWxN8\nrmpS73QzO2Bms5ppt8X6IiKpRrPLEyAUGs1TT+VG9nft2sXy5cs5duwYNTU1LFu2jMsuu4x+/Zo+\nuykiSaAWmOXu68zsVGC1meUHxxa6+8Jm6v0YeKmVtluqLyKSUtSTmQSGDh3KpZdeSlpa+F9HVVUV\nS5cupaqqKsGRiUhT7r7T3dcF21XAZmBkcNhi1TGz64CtwKZWmo9ZX0QkFSnJTBJZWVl86lOfiiSa\nNTU1PPXUM/zzP3+fyZNztWi7SBIyszHABGBlUHSrma0zs1+b2cDgnH7Ad4A8Wk8io+tndFDYIiKd\notUk08weMbMyM1sfVbbAzDYHN8P/NrMBUcdmm9l7wfErOirwrmjkyJFceumlpKenU16+i7lzV7N4\n8WwKCvJ4+ulvM2XKz5RoiiSJYKj8OeD2oEfzIWCsu08AdgL3B6fOBR5w95r6qs002bS+hs1FJKW1\n5ZnM3wA/A56IKlsC3O3udWb2Q2A2MNvMzgGmAuOBbGCpmZ3pekl3m5122mnk5OTwpS/dTVnZzzh+\n0fb7Gz3PKSKdz8x6EE4wn3T35wHcfVfUKQ8DLwbbHwe+ZGYLgEzgmJkddPeHottsoX4jRUVFzJ07\nN7Kfk5NDTk7OSV2PiHQ/BQUFFBQURPaLi4vj/h2tJpnu/lczG92kbGnU7grgS8H254BF7l4LFJvZ\ne8BEGoaSpA2GDRuG+3BiLdpeWlqXiJBEpLFHgUJ3f7C+wMyy3H1nsPtFYCOAu18adU4ucKBpgtlS\n/aZCoRC5ufqhKSInp+kP1Ly8vLh/RzyeybyRhhmTI4HtUcdKaHggXk7AmDG9ibVoe58++1DHsEji\nmNklwA3AZWa2Nmq5oQVmtt7M1gGTgG+1oa2HzezCYPeE64uIJLOTWsLIzL4HHHX3Z+MUjwTmzZvB\nihW5bNmSR7hHM7xo+xVX/AN//vOf+cQnPkHPnj0THKVI9+PurwPpMQ690oa6eU32vxm1Pf3koxMR\nSR7tTjLNbAZwNXBZVHEJMCpqPzsoO46eK2pZ/aLtc+bcT0lJLenpu7n22n9g2LChlJaWsmTJEj79\n6U+TkaEJqNK9dcZzRSIicuLammQaUTMig6Ghu4BL3f1w1HkvAE+b2QOEh8nPAFbFalDPFbUuetF2\nd+ftt99m8+bNABw4cIAlS5YwceJERo8e3VIzIl1aZzxXJCIiJ67VJNPMngFygMFm9gGQC9wD9ATy\nzQxghbvf4u6FZrYYKASOArdoZnl8mBkTJkxg0KBBrFy5ktraWmpra/nDH57npZfe4/DhgWRn92De\nvBmEQko6RUREJLHaMrv8+hjFv2nh/PnA/JMJSpp3+umnk5GRwV/+8he2bNnK/PmbKCtbSP1zmytW\n5JKfP1OJpoiIiCSU3viTgjIyMrjyyit56aV3oxJMqF9L8/vfb/Y3gIiIiEinUJKZok455RSOHBlM\nrLU0N27cRUVFRSLCEhEREQGUZKa0kSPTiLWW5qmn1pCfn8/q1as5evRoIkITERGRbu6k1smUxIq9\nluadTJ16AQDvvvsuH3zwAYMGDeFXv3qN0lJn5Mg0TQ4SERGRDqckM4VFr6VZWlrHiBFpzJ59O7t2\nlbFjxw4APvhgO//6r69ocpCIiIh0KiWZKS56Lc0G49m2bRtr165l8eK1lJX9nKaTg+bMuT9GPRER\nEZH40DOZXdTo0aO55pprOHJkELEmB23evJeqqqpEhCYiIiLdgJLMLqxHjx6cffYAYk0O6t17H3/8\n4x9ZuXIl+/fvT0R4IiIi0oVpuLyLiz05aBZTp15AXV0dW7duZevWrWRnZ9O376k88MALlJTUaYKQ\niIiInBQlmV1crMlBd9xxM/v27aW8vDxy3po1a/X2IBEREYkbJZndQOzJQVBRUUFhYSElJSXNThC6\n++75/Nd/3dep8YqIiEjq0zOZ3diQIUO49NJLufrqqzl0KIPm3h702muvsXXrVo4ePUpR0TamTctj\n8uRcpk3Lo6hoWyJCFxERkSSnnkwhIyOD8eMzWbOmmsaJZjUDBx6hvLyc8vJy/vjHl/nRjwopKVmA\nhtRFRESkJerJFCA8QWjcuFwaZqJXM2LEd5g69cLIOYsWvRWVYEL9kPr3vvdo5wYrIiIiSU89mQLE\nniA0b953yMoaxrZt2ygqKqKyshexhtTXry/jT3/6EyNGjCArK4uKij3MmfOYZqmLiIh0Y0oyJaK5\nCUIf+chH+MhHPsKjj/6VwsLjh9QzMg5TWlpKaWkp5eW7+OEPN7Nz5/1oSF1ERKT70nC5tNmCBTcd\nN6R+2mnfZurUCyLnLF68NirBhPoh9Vtu+THbtm2jpqYmcq4mEYmIiHRd6smUNos9pH43Q4cOZseO\nHezYsYN9+3oTa0j9gw8O87e//Q2Avn37cuRILXfe+Se2b5+PejxFRES6HiWZckKaG1I/88wzOfPM\nMzn//GVs2hR7lnq9mpoafv7zJWzffvy6nDNnzuGXv7yTjIwM+vXrR3HxB3q+U0REJAUpyZS4uu++\nf2HlysavsRwz5nvMmfMVevU6hYqKCmpra5udRFRUVM2f//xnACoq9jB//kZ27Gh4vvNvf5vDkiW3\ncsYZYyO1ioq2KREVERFJMkoyJa5iD6l/K5L0uTuVlZU8++zqmJOIons8Fy16ix07Gvd2FhXN42tf\nu5W77rqOU089lf37q7jjjtfYvv2HNAy7/xv5+bc1SjSViIqIiHQuJZkSd80NqQOYGZmZmTz44EwK\nCxv3eJ5++my+850vcuqpfdm/f3+zvZ2VlT05ePAgBw8ebGbY/V5mzLid731vKn379mXv3n383//7\naqPnP994Yw75+bcxduyYRq0rGRUREYmPVpNMM3sEuAYoc/fzgrJM4L+A0UAxMNXd9wXHZgM3ArXA\n7e6+pGNCl1QWu8fzzkYJ3eLFa2P2dmZmNvR2NpeIVlSks3v3bnbv3h0zEd26dR7Tp9/KnXd+jl69\netG7d2/27KnkW99azocfNvSK/vWv3+O///trnHnmOHr27EnPnj3Ztm17q4loW5NVJbWpx8yygSeA\n4UAd8Ct3/5mZ5QLfBMqDU+9x91ei6p0ObAJy3X1hjHabva+KiKSitvRk/gb4GeGbar27gaXuvsDM\nvgvMBu42s3OAqcB4IBtYamZnurs3bbSoqOikg0+UgoICcnJyEh3GCUu2uFvq8QT4f//v67z5Zn1v\n55vAxxg3LpfHH/83hg8fSlVVFc8++1arw+7NJaJ79/bk0KFDHDp0iH379vHzny/hww8bJ6Pbtv2A\n2267lVtvvQKA8vJdzJ9fSFnZj6lPRJcvv4sHHriM00/PpkePHpSVlfN//s/LQbIajvuvf/0ev/3t\nVxk3LkR6ejrp6els317CNdf8iq1b76WlGfZtSUTjmdTWn7Nx41bOPXfsSX1fInTCvaUWmOXu68zs\nVGC1meUHxxbGSiADPwZeaqHdmPfVpiel8r2zLZLtPhVvur7U1ZWvDTrm3tJqkunufzWzpv/3uA6Y\nFGw/DhQQvhl+Dljk7rVAsZm9B0wEVjZtt7i4uP1RJ1iq/kVLtbijeztff30Zl1xyGfPmNSRgp556\nKg8+eNtxw+5jx/4bv/zlXQwfPpSamhpCoTdbTUSh+WS0srJnZG/x4rWUlTVOREtL/50HH2xIRBsn\nqwVADtu2/YA77mg4p/68rVuPn2H/ta/dzv9v72xjpCqvOP77D75haxWxLArKgoZoExM0WEzQSmo1\ntB+gaWJiNLZKTA0ttkHRIk1ZbH0BrRBrQz+gJII1gtagbZUiFmIaRbS8CIKKcZfldUWradWAAqcf\nnmdhZudlZ4e5O/duzi+ZzJ1nnnv3f5575+y5z8u506ZNoF+/fuzb9xEzZrzOnj0PHrFv9eo7efjh\nKxk6dAi5XI7du/cyZcqKgh7YV165i4ULJ3L22UPJ5XLkcjl27drDDTc8TXv7feT31C5dej3Dhp2D\nJNrbd3LNNYtpa7sHeJCNG6fx6qu/Ydmym2huHoYktm/fwYQJj9La+jvypx+88MJkRoxoJpfLIYnW\n1u3MnPk4u3db3QLfaoLapH2Lme0F9sbtzyRtBYbEr1VqH0kTgQ84mmS2FOX8agFZ9p3VkDU/1VPc\nvuzSl22DZHxLrXMyB5lZBwSHK2lQLB8CvJZXbxdHna/j9JjO3s5Zs4xZs4p7PUsPuxcu+nnkkV/y\nzjuFgei557awaNFMzjprMPv37+fAgQMsWbKuZDDa1CROOeUUvvzyy6oC0WrqVKr38cf9+PzzEIss\nWLC6aPHTrl0PMHduuaA21Nmx435aWoqD2vb24p7aqVMLj9XWVrzYavLkwjqtrcXTDyZN6trj+zYd\nHRhGZusAAAhZSURBVHOPtPnKlbcxY8aFNDV9Ewjzczs69nHvvZsKeoZXrrydlpZRNDUNinU+ZNas\n9QV1Xn75du6+ezRnntl05Fj5772FpGZgFOFG+jJgiqQbgDeBaWb2qaSvAXcCVwF3VDhcOb/qOI6T\nSeq18KdoONxxeovuht1LB6JHe0T79+8PwLx5U9i8uTgYfeyxu47UffbZTSUD0ZEjT2XMmDEcPHiw\nbM/p4ME5BgwYwKFDhzh06BADBx4kdGz1fKi/nkFtT49Va49vR8dcliwpDHyfeurNEvUe4oknjtZb\nvHhNUZ29ex/i8ccLj9XbxKHyZwhzzz+TNB/4rZmZpHuA3wM3A7OAeWb2RQyCq42E3a86jpNpVGK6\nZHGlMFz+17yFP1uBcWbWIWkwsMrMLpA0HTAzmxPrLSdMci8aLh8/frzt37//yOfm5maGDx9eF6OS\nprW1NTNa88mqbsiu9qzqhuxob21tLRjmOemkk1i+fHmiXZqSjgP+BrxoZg+X+P6Iz5T0CmGOOsAA\n4BAw08zmd9mnpF/teuws+85qyMp1VytuX3bpa7b1hu+sNshsJjjMC+PnOcB/zGxOnKA+wMw6F/78\nGRhDGCZ/CSi58MdxHCerSFoEfGRmt+WVDY7zNZE0FbjEzK7rsl8L8L8yq8tL+tVEDXEcx0mQalIY\nPQmMAwZKagdagNnA05ImAdsJK8oxsy2SlgJbgK+An3mA6ThOX0LSWOB6YJOk9YRh7RnAdZJGEdIa\ntQG3VHGsBcCfzGwdMAdY2tWvOo7jZJWqejIdx3Ecx3EcpyfkGvFHJY2X9I6k9+KwUGqR1CZpo6T1\nktbGsgGSVkh6V9I/JJ3aaJ0QEudL6pD0Vl5ZWa2S7pK0TdJWSQ1bQVFGd4uknZLWxdf4vO9SoTtq\nGSrpn5LelrRJ0i9iearbvYTuW2N5qttd0omSXo+/x7cl3RfLU93e5Sh17Zeo84eof0PsKe0sL+lH\n0+SfErLvgXguN0j6i6RvJG1HGd11ty3v+9slHZZ0elL6uyMp+yTdGs/fJkmzk7ShEgldm5dIWhv9\n01pJo5O2oxw12HdRd/vW5FvMrFdfhMD2fcJTLY4HNgDn97aOHuj9gDA3Kr9sDnBn3P4VMLvROqOW\nywjpVN7qTivwLWA9YcpEczwnSpHuFkLC6651L0iL7qhnMDAqbn8deBc4P+3tXkF36tsdODm+9wPW\nAGPT3t4VbCm69rt8/33g73F7DLAmbpf1o2nyTwnZ9z0gF7dnA/f3Fdvi90OB5UArcHofO3fjgBXA\ncfHzGX3MvlXA1Xn7r8qafZX2rcW3NKIn89vANjPbbmZfAU8RkhCnFVHc4zuRkCyZ+P7DXlVUBjP7\nF/BJl+JyWo8kzjezNqAzcX6vU0Y3lE71MpGU6IaQz9DMNsTtz4CthH8SqW73MrorJRRPTbub2Rdx\n80TCb/MTUt7e5ahw7Xcykfi0NQtZOk6V1ERlP5oa/5SEfWa20swOx/3XcHTlfq+S0LkDmEflfKq9\nQkL2TSYEJgfjfh8lpb87ErJvD9DZu3caIVd4QzgG+yrt22Pf0oggcwiwI+/zTtKdsN2AlyS9Ienm\nWNZkeUmTgTQnTR5URmvX85DGxPlTYjf+o3nd8qnVraOJuddQ/hpJnf483Z2pxlLd7pJyCgtu9gKr\nzWwLGWrvHlLOX1byo1nyT7XYl88k4MXE1B0bPbZN0gRgh5lt6i2Rx0At524k8B1JayStauRwchXU\nYt90YK7CIukHCI+GTSu1+MZy8URZGjInM2OMNbOLgR8AP5d0OcVJkrO0eiorWucDI8xsFCGYeKjB\neiqiLom5ycg1UkJ36tvdzA6b2UWEHqzLJY0jI+1dB2rJYZeltqjaPkm/Br4ysycT1FNPKtomqT8h\nS0H+kyV69xFWx0Y1Wo8jTD+7lPAUrKXJSqor1dj3GHCrmZ0DTAUWJiup4XTrWxoRZO4Czsn7PJQG\ndil3h5ntie/7gGWErvKOzm5lhaTJHzZOYbeU07oLODuvXqrOg5ntszjxA1jA0SHO1OlWSMz9DLDY\nzJ6Lxalv91K6s9TuZvZf4AVgNBlo7xopp7+SH92bIf9Ui31IupFw41+QhzRl9NS2cwnzhjdKao3l\n/1Z6Hy9ay7nbCTwLYGZvAIclDUxeak3UYt8YM1sGYGbPkKKpOSWoxTf2OPZpRJD5BnCepGGSTgCu\nBZ5vgI5ukXRy7OlB4fnDVwObCHpvjNV+AjxX8gCNQRTecZXT+jxwraQTJA0HzgPW9pbIEhTojhdw\nJz8CNsfttOmGcLe6xQqf/JKFdi/SnfZ2l3RG5xB+7Pm5irCwJwvtXY6uv9l8ngd+DCDpUuDTOFxV\nyY+mzT/V1T6FjAd3ABPM7EDS4ruhbraZ2WYzG2xmI8xsOCEgu8jMGnmTUO9rcxnw3bjPSOB4M/s4\nQf3dUS/7On9j2yRdEfe5EngvSfFVUIt9lfbtuW/pbmVQEi9gPGE16zZgeiM0VKlzOGHl2HpCcDk9\nlp8OrIw2rABOa7TWqOtJYDdwAGgHbiI8xq6kVsJ8kfcJiz6uTpnuRcBbsf2XEeaZpUp31DKW8JjA\nzutkXby+y14jadBfQXeq2x24MGpdD2wEpsXyVLd3BXtKXfu3AD/Nq/PHqH8jcHFeeUk/mib/lJB9\n2wjJ6tfF1/y+YluX439AY1eXJ3HujgcWE/6fvglc0cfsG02Y274eeI1wk5BF+4r2jeU99i2ejN1x\nHMdxHMepO77wx3Ecx3Ecx6k7HmQ6juM4juM4dceDTMdxHMdxHKfueJDpOI7jOI7j1B0PMh3HcRzH\ncZy640Gm4ziO4ziOU3c8yHQcx3Ecx3HqjgeZjuM4juM4Tt35PwVnQUU4Pt85AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106bf0e80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tt = sc.arange(0,p[\"timeMax\"],p[\"timeStep\"]*p[\"timeMax\"]/3.0)\n",
    "fig=gr.figure(figsize=(11,5))\n",
    "rows=1; cols=2\n",
    "ax1=fig.add_subplot(rows,cols,1)\n",
    "ax2=fig.add_subplot(rows,cols,2)\n",
    "ax1.plot(p[\"sampTimes\"], uSol[0],'k',lw=3, alpha=0.35, label=r'($t,y(t)$)')\n",
    "ax1.plot(tt,f(tt,p),'bo',label=r'($t,u(t)$)')\n",
    "ax2.plot(p[\"sampTimes\"], uSol[0],'k', lw=3, alpha=0.35, label=r'($t,y(t)$)')\n",
    "ax2.plot(tt,f(tt,p),'b',label=r'($t,u(t)$)')\n",
    "ax2.set_xlim(1,1.01); ax2.set_ylim(254,258)\n",
    "ax1.legend(); ax2.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}