{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "ZyYxarWXICGB" }, "source": [ "\"Licencia
Python en ciencias e ingeniería: tutoriales basados en ejemplos por Sergio Gutiérrez Rodrigo y Adrián Navas Montilla se distribuye bajo una Licencia Creative Commons Atribución-NoComercial-CompartirIgual 4.0 Internacional." ] }, { "cell_type": "markdown", "metadata": { "id": "cU67irCkAoSz" }, "source": [ "---\n", "Programa de Recursos en Abierto en la UZ (2022)\n", "\n", "**Python en ciencias e ingeniería: tutoriales basados en ejemplos**\n", "---\n", "Universidad de Zaragoza\n", "---\n", "*PRAUZ-739*\n" ] }, { "cell_type": "markdown", "metadata": { "id": "aZ4I0EbyAqVW" }, "source": [ "##
Resolución numérica de la ecuación del calor
\n", "\n", "\n", "-----------------------------------------" ] }, { "cell_type": "markdown", "metadata": { "id": "ctgg3tv4d4a9" }, "source": [ "## Introducción al método de las diferencias finitas\n", "\n", "\n", "### Derivadas discretas (aproximadas) \n", "\n", "La aproximación por diferencias finitas es el más antiguo de los métodos aplicados para obtener soluciones numéricas de ecuaciones diferenciales. La primera utilización de este método se atribuye a Leonhard Euler en 1768. La idea de los métodos de diferencias finitas es en realidad bastante simple, ya que corresponde a una estimación de una derivada mediante la diferencia discreta del valor de la función en dos puntos cercanos.\n", "\n", "La derivada de una función $u(x)$ en $x$ se define como:\n", "\n", "$$\n", "\t\\frac{\\partial u }{\\partial x} = \\lim_{\\Delta x \\rightarrow 0} \\frac{u(x+\\Delta x) - u(x)}{\\Delta x} \\,.\n", "$$\n", "\n", "\n", "Si quitamos el límite de $\\Delta x \\rightarrow 0$ en la expresión anterior, considerando un $\\Delta x$ finito, entonces obtendremos el concepto de derivada discreta como una diferencia finita. Ésta se calculará como el ratio de la diferencia de la función evaluada en $x+\\Delta x$ y $x$ respecto de $\\Delta x$. \n", "\n", "En el caso de utilizar una diferencia finita, se cometerá cierto error y éste se puede derivar utilizando una expansión en serie de Taylor:\n", "\n", "$$\n", "u(x+\\Delta x )=u(x)+\\left(\\frac{\\partial u}{\\partial x}\\right)\\Delta x+\\left(\\frac{\\partial^2 u}{\\partial x^2}\\right)\\frac{\\Delta x^2}{2} +\\left(\\frac{\\partial^3 u}{\\partial x^3}\\right)\\frac{\\Delta x^3}{6} + ...\\,\\,,\n", "$$\n", "\n", "Despejando la primera derivada de la expresión anterior obtenemos:\n", "\n", "\n", "$$\n", "\\left(\\frac{\\partial u}{\\partial x}\\right)= \\frac{u(x+\\Delta x) - u(x)}{\\Delta x} -\\underbrace{\\left(\\frac{\\partial^2 u}{\\partial x^2}\\right)\\frac{\\Delta x}{2} -\\left(\\frac{\\partial^3 u}{\\partial x^3}\\right)\\frac{\\Delta x^2}{6} - ...}_{\\mbox{error de truncamiento}}\\,\\,,\n", "$$\n", "\n", "lo que permite expresar la aproximación de la primera derivada como:\n", "\n", "$$\n", "\t\\frac{\\partial u }{\\partial x} \\approx \\frac{u(x+\\Delta x) - u(x)}{\\Delta x} \\, .\n", "$$\n", "\n", "En 1D, vamos a considerar un dominio espacial $[0,L]$ que se discretizará en $N=L/\\Delta x+1$ nodos localizados en las coordenadas $x_i=i\\Delta x$, con $\\Delta x$ el paso de discretización espacial. En los puntos de la malla, se aproximarán las variables como $u_i \\approx u(x_i)$\n", "\n", "\n", "Generalmente, se usan dos tipos de aproximación de las derivadas espaciales en función de qué información cojamos:\n", "\n", "- Diferencias regresivas (hacia atrás):\n", "\n", "$$\n", "\t\\frac{\\partial u }{\\partial x} \\approx \\frac{u_{i} - u_{i-1}}{\\Delta x} \\, .\n", "$$\n", "\n", "- Diferencias progresivas (hacia adelante):\n", "\n", "\n", "$$\n", "\t\\frac{\\partial u }{\\partial x} \\approx \\frac{u_{i+1} - u_{i}}{\\Delta x} \\, .\n", "$$\n", "\n", "En el caso de tener que aproximar una derivada segunda, podemos combinar ambas para obtener:\n", "\n", "$$\n", "\t\\frac{\\partial ^2 u }{\\partial x^2} \\approx \\frac{\\frac{u_{i+1} - u_{i}}{\\Delta x} - \\frac{u_{i} - u_{i-1}}{\\Delta x}}{\\Delta x} = \\frac{u_{i+1} -2u_{i} +u_{i-1}}{\\Delta x^2} \\, .\n", "$$\n", "\n", "\n", "### Resolución numérica de la ecuación del calor\n", "\n", "Ahora vamos a utilizar las expresiones anteriores de derivadas discretas para resolver la ecuación del calor en 1D sin generación:\n", "\n", "$$\n", "\\rho c\t\\frac{\\partial T }{\\partial t} = \\frac{\\partial }{\\partial x}\t\\left( k\\frac{\\partial T }{\\partial x} \\right) \\, .\n", "$$\n", "\n", "Asumiendo que $\\rho$, $c$ y $k$ son constantes, podemos reescribir la ecuación de forma compacta como\n", "\n", "$$\n", "\t\\frac{\\partial T }{\\partial t} = b\t\\frac{\\partial ^2 T }{\\partial x^2}\\, .\n", "$$\n", "\n", "con $b=k/(\\rho c)$ la difusividad térmica. El problema se define en un dominio $[0,L]$ hasta un tiempo $t$, considerando también una condición inicial $T(x,0)=T_0(x)$ y unas condiciones de contorno $T(0,t)=T_L$ y $T(L,t)=T_R$. \n", "\n", "Para resolver el problema, el primer paso es expresar la ecuación diferencial de manera discreta, utilizando una la aproximación progresiva para el tiempo y la derivada central para el espacio:\n", "\n", "$$\n", "\\frac{T_{i}^{n+1}-T_{i}^{n}}{\\Delta t}=b\\frac{T_{i+1}^{n}-2T_{i}^{n}+T_{i-1}^{n}}{\\Delta x^2}\n", "$$\n", "\n", "donde $n$ y $n+1$ son dos instantes consecutivos de tiempo tal que $t^n=n\\Delta t $, mientras que, $i-1$ y $i$ son dos puntos vecinos de la coordenada $x$ discretizada. Si se dan las condiciones iniciales, entonces la única incógnita en esta discretización es $u_i^{n+1}$. Por lo tanto, podemos actualizar el valor de la temperatura en cada punto utilizando la siguiente expresión:\n", "\n", "$$\n", "T_{i}^{n+1}=T_{i}^{n}+ b\\frac{\\Delta t}{\\Delta x^2} \\left(T_{i+1}^{n}-2T_{i}^{n}+T_{i-1}^{n}\\right)\n", "$$\n", "\n", "Pero... \n", "\n", "- ¿Con qué valores de $T_{i}$ empezamos en el tiempo $t=0$, es decir, en $n=0$? Impondremos la condición inicial:\n", "$$T_{i}^{0}=T(x_i)$$\n", "\n", "- ¿Qué ocurre cuando queremos caluclar la temperatura en el primer nodo ($i=0$) o en el último? Entonces no podemos utilizar la expresión general anterior ya que nos falta información de lo que hay al otro lado del contorno. En este caso, impondremos las condiciones de contorno. En el caso más sencillo, cuando nos dan la temperatura en los extremos, tendremos condiciones de contorno de tipo Dirichlet y bastará con hacer:\n", "$$T_{0}^{n+1}=T_L$$\n", "$$T_{N-1}^{n+1}=T_R$$\n", "donde $T_L$ y $T_R$ son los valores de temperatura en los extremos izquierdo (L, left)y derecho (R, right), que son dato del problema.\n", "\n", "Y una última cuestión... debemos elegir el valor de $\\Delta t$. Para que el método numérico funcione adecuadamente y no se vuelva inestable, este valor no se puede escoger de forma arbitraria. Lo calcularemos de la siguiente manera:\n", "\n", "$$ \\Delta t=\\sigma \\frac{\\Delta x^2}{b}$$\n", "\n", "donde $\\sigma$ es una constante que debe estar dentro de unos límites. A esto lo denominamos *condición de Courant–Friedrichs–Lewy (CFL)*." ] }, { "cell_type": "markdown", "metadata": { "id": "LITVC6KmPNzZ" }, "source": [ "## Trabajo práctico\n", "### Primer paso: cargar librerías de Python\n", "\n", "Para comenzar, debemos cargar todos los paquetes de Python necesarios: " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "FNuY3zZTdbtN" }, "outputs": [], "source": [ "import numpy as np # Librería para cálculo numérico\n", "import math # Librería para utilizar símbolos matemáticos como el número pi, que se escribe como math.pi\n", "from scipy import special\n", "import matplotlib.pyplot as plt # Librería para poder dibujar gráficas\n", "import matplotlib.animation as animation\n", "from matplotlib.animation import FuncAnimation\n", "from IPython.display import HTML" ] }, { "cell_type": "markdown", "metadata": { "id": "xzTBZ3ZZQGR3" }, "source": [ "### Apartado 1: Primeros pasos para la resolución de la ecuación del calor\n", "\n", "Para comenzar, se pide resolver el siguiente problema\n", "\n", "$$\n", "\t\\frac{\\partial T }{\\partial t} = b\t\\frac{\\partial ^2 T }{\\partial x^2}\\, .\n", "$$\n", "\n", "con $b=0.3$ m$^2$/s en un dominio $[0,2]$ m hasta un tiempo $t=0.2$ s, dada una condición inicial \n", "\n", "$$ T(x,0)=\\left\\{ \\begin{array}{l}\n", "2 \\hbox{ ºC si } x<1 \\\\\n", "1 \\hbox{ ºC si } x\\geq 1 \n", "\\end{array}\\right.$$\n", "\n", "y unas condiciones de contorno $T(0,t)=2$ ºC y $T(L,t)=1$ ºC. Para ello, se utilizará el esquema en diferencias finitas expuesto arriba, que nos permite calcular la evolución de la temperatura a o largo del tiempo.\n", "\n", "Se pide lo siguiente:\n", "\n", "1.a) Completa el código inferior para resolver el problema propuesto, utilizando las expresiones anteriores.\n", "\n", "1.b) Completa la gráfica y observa el resultado.\n", "\n", "1.c) Prueba distintos valores de $N$ y de $b$ y observa la solución. \n", "\n", "1.d) Prueba $\\sigma=0.05$ y $\\sigma=0.9$ y explica qué ocurre. Encuentra (probando) el valor aproximado de $\\sigma$ para el cual la solución deja de ser estable y se vuelve oscilatoria.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "lV7vNgQN_hII", "outputId": "3ef2385f-6de2-452c-da30-1a3dedcb85ed" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "La longitud del dominio es: 2.0\n", "El numero de nodos espaciales es: 20\n", "El numero de pasos temporales es: 18\n", "El paso de tiempo es es: 0.011080332409972297\n", "El tiempo final es: 0.18836565096952906\n", "\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAE9CAYAAAAGZmUpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7RElEQVR4nO3dd3hUVfrA8e/JZJKZtAkQCCWhiHRIQIIiva6IgGUp1gXrz7aK7qqsuoisq9gWde1t7aK4YkFlpQqKCAHpCAIihBJ6ekLK+f1xQ0iZJJNk7tzM5P08zzzJ3HPLy3idN+eeprTWCCGEaLiCrA5ACCGEtSQRCCFEAyeJQAghGjhJBEII0cBJIhBCiAZOEoEQQjRwwVYHUFMxMTG6bdu2VochhBB+Ze3atUe11k3dlfldImjbti3JyclWhyGEEH5FKfV7ZWXyaEgIIRo4SQRCCNHASSIQQogGzu/aCIQINPn5+aSkpJCbm2t1KCIAOBwO4uLisNvtHh8jiUAIi6WkpBAZGUnbtm1RSlkdjvBjWmuOHTtGSkoK7dq18/g4eTQkhMVyc3Np0qSJJAFRZ0opmjRpUuPapSQCIeoBSQLCW2pzL5mWCJRS8UqppUqpbUqpLUqpO93so5RSzymldiqlNiqlzjElmI0fw+zuMCPa+LnxY1MuI4Q/OnToEJdffjnt27ena9eujB49mh07dtT5vMuWLWPMmDEAfPHFF8yaNavK/fv161fnawJMnz6dRYsWVbmPt+KJiIioUWz1lZltBAXAX7TW65RSkcBapdRCrfXWUvtcCHQofp0HvFT803s2fgxf3gH5Ocb7tH3Ge4CEiZ6fY/FMSEsBVxwMn+75sULUY1prLr30UiZPnsycOXMAWL9+PampqXTs2NFr1xk3bhzjxo2rcp+VK1d65VozZ86sV/H4A9MSgdb6IHCw+PcMpdQ2oBVQOhFcDLyjjWXSVimlopVSLYqP9Y7FM88kgdPyc8iZfx/LD4Wjg8MZdU57CAlnVUoev53ML7Nr+0Nfc+6mGWUSScFnf2b1rmP83uoiwkODGZfYEoClvxzmUHrZZ3PRTjsX6hWweCY6LYVsR3M2dLqT31tdBEBMRCgju8YC8NXGg6Tnlr1+c5eDoZ2aAfD5+v1knyosUx7fKIwBHWIA+GRtCvmFRWXK28WE0/esJgDMWb2X8uvRdYyNoHebxhQUFjF3bUqFj69riygS46PJzS9k3s/7K5QnxLno1tJFRm4+8zdW/M/Wu00jOsZGciLrFAu2HKpQfl67xpzVNILDGbks3na4QvmAs2OIbxzGgZM5fLfjSIXyIZ2a0sLlZO+xbH7YdbRC+YgusTSNDGXXkUxW/3a8Qvno7i1whdnZfiiDdXtPVCgfl9iS8NBgNu9PY0dqBqN7tMBht1XYz18tXboUu93OzTffXLKtZ8+egJEk7r33Xr755huUUjz44INMmjSJZcuWMWPGDGJiYti8eTO9e/fmvffeQynFggULmDp1KjExMZxzzpkK/ltvvUVycjLPP/88qamp3HzzzezevRuAl156iX79+hEREUFmZmatrlvalClTGDNmDOPHj6dt27ZMnjyZL7/8kvz8fObOnUvnzp1rFE9mZiYXX3wxJ06cID8/n0ceeYSLL77Y5P8yvuWTXkNKqbZAL+CnckWtgH2l3qcUbyvzjaKUugm4CaB169Y1u3haxS83AOep41yw8krjzXLjR1/gHG0jh1CycJCtHbQOOoxRuTkjuCiXLusf4T9rDmOLbMq4+FEQ3pQ3Vuzm+13Hyux7c6O1XFjwIuTnoIDw3IP0XD+dOWv28kXRAHq3aVSSCJ5ZtINfD2eWOX5Qx6YlieDxb37hQFrZRDO6R/OSRDDzyy2k55aNdULvuJJE8OBnmykoKpsKpvRraySCIs3fPt1U4XO6bWh7EuOjycorcFt+76hOdGvp4kRWvtvymRd3o2NsJAfTct2WPz0hkbOaRrD3WLbb8pev7k184zB2pGa4LX/3+nNp4XKyaX+a2/JOt0bSNDKUtXtOuC3v3aYRrjA7P+w8ysz5WyuUD+rYlPDQYJZtP8xT3+4gJDiIMQktK+znr05/obrz6aefsn79ejZs2MDRo0fp06cPgwYNAuDnn39my5YttGzZkv79+/PDDz+QlJTEjTfeyJIlSzj77LOZNGmS2/PecccdDB48mHnz5lFYWEhmZtl7vqbXHTBgQJX/xpiYGNatW8eLL77IU089xeuvv16jeBwOB/PmzSMqKoqjR4/St29fxo0bF1DtOqYnAqVUBPBfYKrWOr18sZtDKiyirLV+FXgVICkpqWaLLLvijMdB5RSGNSV95GxUfhbRwafgVBY5WekU5mYSnJ9FdH4WjU5lY//1C7enbaQyeS3kX5AHPPc3AN4NdlDUNIbCsBiKil+OX7+qUCMJU6f4V5MvuP+Gh7DbznwEH9zYl8JyX9QhwWeacT67vT9FZf/gJ7RU+cK7B1N+CWpnqb9ev79vWIV/hzPEVnKeVX8bXqE8LNQojw4LcVse4TBuoZbRDrflkcXlHWIj3Ja7nEZf5x5xLrfl0WFGed+zmlRZPrxLM7fljcNDABiT2IJBHSvOt9Ukwiif1Cee0T1aVCiPKS6/uGcrnvp2B8ezTlXYx9smvfJjhW1jElpwzfltyTlVyJT/rK5QPr53HBOS4jmedYpb3ltbpuyj/zu/VnF8//33XHHFFdhsNmJjYxk8eDBr1qwhKiqKc889l7i4OMCoQezZs4eIiAjatWtHhw4dALj66qt59dVXK5x3yZIlvPPOOwDYbDZcLledrltdIrjssssA6N27N59++mmN49Fac//997N8+XKCgoLYv38/qampNG/evNrP0F+YmgiUUnaMJPC+1rrifwGjBhBf6n0ccMCrQQyfXraNAMDuxDbqURoljC2zq9Pd8bO7u00kRLaAK+ZA1lHIOgJZR1BZR7BlHcVW/J5j2+BUhtuwgjNSaP7DdIjtCs26QbMuNI2spOGpuI2iWTVtFLFRjko+BENzV+XlSqkqy21BVZcH24KqLLdXUx4abKO5q/JHLg573crDQoIJC6n8dg8PDSY8tPLyZlGhAKTn5Fe6jz/q1q0bn3zyidsyXf6vilJCQ0NLfrfZbBQUGDVRb/yVXJvrVuX0MZ7uX97777/PkSNHWLt2LXa7nbZt2wbc4D/TEoEy7og3gG1a639VstsXwO1KqTkYjcRpXm0fgDNfmLVt7K0kkTByJrTsWf3xlSUSWwisfx9OlaqGNmprJIXYbmcSxP518NXUujV2izoLDbbhsAeR5oNEUNVf8M4QW5XljcNDalQDGDZsGPfffz+vvfYaN954IwBr1qwhOzubQYMG8corrzB58mSOHz/O8uXLefLJJ/nll1/cnqtz58789ttv7Nq1i/bt2/Phhx+63W/48OG89NJLTJ06lcLCQrKysoiKiiopr+l166q6eNLS0mjWrBl2u52lS5fy+++VTuLpt8ysEfQHrgE2KaXWF2+7H2gNoLV+GfgaGA3sBLKBa02JJGFi7b80zUokY5+D7uMhbS+kboHUrXB4i/H7jm9AF1V+zvwcIx5JBD6V0Cq65FFWoFBKMW/ePKZOncqsWbNwOBy0bduWZ555hkGDBvHjjz+SmJiIUoonnniC5s2bV/qF7HA4ePXVV7nooouIiYlhwIABbN68ucJ+zz77LDfddBNvvPEGNpuNl156ifPPP5O8Lr300hpdt66qi+eqq65i7NixJCUl0bNnTzp37mxKHFZSVVXD6qOkpCTtd+sR1LT7aX4uHN1uJIXPbql8vzs3GLUI4de2bdtGly5drA5DBBB395RSaq3WOsnd/jLXkC/UtEZid0CLROO19FH3j5YAnk2EludAt0uh2yUQXcMeVUIIgUwxUf8Nn248SirN7oQLHoURDxuPkBb+HZ7pAa+PgB9fqLTLrKibxxf8wu0frLM6DCG8TmoE9V11bRQDpsLx3bDlM9gyD/53v/GKPw+6XQZdL4Y9K2RktBccSstl/b6TVochhNdJIvAH1T1aanwWDLzbeB3dCVvnGYlhwX3GSwWdaXyWXke1FuUI9kmvISF8TR4NBZqYs2HQPXDLD3DbGgiNqtgD6XSvI1EjLqedzLwCior8q4OFENWRRBDImnaEPPcD2kjbB4Xy121NRDntaA0ZuTUflCQq99FHH7Fnzx6rw2jQJBEEOldc5WUvD4Cdi30Xi59r0ySc89o15lRhFWM8/NA///lPunXrRkJCAj179uSnn8pPCVbWkCFDqG0X7vJTO7/33nvs3buXtm3b1up8tXHy5ElefPHFWh9//PhxRo4cSYcOHRg5ciQnTlScrBDguuuuo1mzZnTv3r3M9vXr19O3b1969uxJUlISq1cbU4bk5+czefJkevToQZcuXXjsscdqHWNNSSIIdJX1Oup7KxTkwXuXwQeXw7Fd1sTnR0Z2jeWj/zufppGh1e/sJ3788Ufmz5/PunXr2LhxI4sWLSI+Pr76A2up/NTOV199Nffcc49p13Onrolg1qxZDB8+nF9//ZXhw4dXuq7BlClTWLBgQYXt9957Lw899BDr169n5syZ3HvvvQDMnTuXvLw8Nm3axNq1a3nllVd8VlOSRBDoEiYao5hd8YAyfo59DkY9Brf9ZEyVsed7eOE8+PZByE2zOmJRHS8utHTw4EFiYmJK5uOJiYmhZUtjdtXFixfTq1cvevTowXXXXUdeXl6F40svzPLJJ58wZcoUAFJTU7n00ktJTEwkMTGxJAGc3l9rzT333EP37t3p0aMHH330EWAsZjNkyBDGjx9P586dueqqq9zOPbRr1y5GjRpF7969GThwYMmo4ylTpnDHHXfQr18/zjrrLLfzKE2bNo1du3bRs2fPWiWhzz//nMmTJwMwefJkPvvsM7f7DRo0iMaNG1fYrpQiPd2YfzMtLa3k81ZKkZWVRUFBATk5OYSEhJSZ6qJ0/F27diUhIYG//vWvNY7fLa21X7169+6thZelH9L6s1u1fsil9RPttU5+S+vCAqujqnd2H8nUQ55cqpdsS/Xqebdu3er5zhs+0vqRWK0fijrzeiTW2F4LGRkZOjExUXfo0EHfcsstetmyZVprrXNycnRcXJzevn271lrra665Rs+ePVtrrfXgwYP1mjVrtNZah4eHl5xr7ty5evLkyVprrSdOnFiyf0FBgT558mSZ/T/55BM9YsQIXVBQoA8dOqTj4+P1gQMH9NKlS3VUVJTet2+fLiws1H379tUrVqyoEPewYcP0jh07tNZar1q1Sg8dOlRrrfXkyZP1+PHjdWFhod6yZYtu3759hWN/++033a1bt5L36enpOjEx0e1ry5YtFY53uVxl3kdHR1f6+Za/ltbGf+/4+HgdFxenW7Zsqffs2aO11vrUqVN60qRJOiYmRoeFhelXXnmlwvmOHTumO3bsqIuKirTWWp84ccLtdd3dU0CyruR7VbqPCoiMhYtfgD43wDfTjO6la16HCx+HNt5ZPjAQhAQH8dvRLA5nmDjz5DfT4FDFdRNKpKyBwnJ/mefnwOe3w9q33R/TvAdc6P7xRUREBGvXrmXFihUsXbqUSZMmMWvWLHr16kW7du1KVimbPHkyL7zwAlOnTvXon2HmVNOZmZmsXLmSCRMmlGwrXVu55JJLCAoKomvXrqSmplYba2RkJOvXr/fo3+UNL730ErNnz+aPf/wjH3/8Mddffz2LFi1i9erV2Gw2Dhw4wIkTJxg4cCAjRozgrLPOKjk2KioKh8PBDTfcwEUXXVSyFGhdSSIQZ7TsBdctgC2fwrfT4T8XGtNXjJwJe1c1+EFpUcVrK1g6lqB8EqhuuwdsNhtDhgxhyJAh9OjRg7fffrtklbLqlJ52uiZTM+s6TDVdVFREdHR0pV/epY+v6jqnZWRkMHDgQLdlH3zwAV27di2zLTY2loMHD9KiRQsOHjxIs2bNqr1GaW+//TbPPvssABMmTOCGG24oudaoUaOw2+00a9aM/v37k5ycXCYRBAcHs3r1ahYvXsycOXN4/vnnWbJkSY2u744kAlGWUtD9j9DxQlj5HHz/DGz90lhCqKj4f8gGOigtIjQYW5AiPcfE7qOV/OVeorJpzV3xcO1XNb7c9u3bCQoKKllMZv369bRp04bOnTuzZ88edu7cydlnn827777L4MGDKxwfGxvLtm3b6NSpE/PmzSMyMhIwd6rpqKgo2rVrx9y5c5kwYQJaazZu3EhiYqJH/+bIyEgyMjLKvK9JjWDcuHG8/fbbTJs2jbfffrvGy1a2bNmS7777jiFDhrBkyZKSz75169YsWbKEq6++muzsbFatWlWhBpaZmUl2djajR4+mb9++nH322TW6dmWksVi4FxIGQ6bBn5MhOORMEjitAQ5KU0pZP7q4sl5gw6fX6nSZmZlMnjy5pPFx69atzJgxA4fDwX/+8x8mTJhAjx49CAoKKrOu8WmzZs1izJgxDBs2jBYtzqzw9uyzz7J06VJ69OhB79692bJlS5njLr30UhISEkhMTGTYsGElU0176v333+eNN94gMTGRbt268fnnn3t8bJMmTejfvz/du3evVWPxtGnTWLhwIR06dGDhwoVMmzYNgAMHDjB69OiS/a644grOP/98tm/fTlxcHG+88QYAr732Gn/5y19ITEzk/vvvL1nF7bbbbiMzM5Pu3bvTp08frr32WhISEspcOyMjgzFjxpCQkMDgwYOZPXt2jeN3R6ahFtWbEY2bFUQBBTNO+jYWi90zdwPdW7mY3K+t185Z42moazqtuWhwZBpq4X2VrPtc5WC1APXkBM8eP5iqLgstCeGGPBoS1XP3OAKgeSL4WY1SCFGRJAJRvQqD0uKgTX/YPt/ottiA5ix68LNNjP3391aHIYRXyaMh4ZnyjyO0hmWz4LtZkJkKE96C0IhKDw8UhUVwMM374wi01mW6YgpRW7Vp95UagagdpWDo32Dss7BrMbw9BjKPWB2V6VxOO+k5+bX6n60yDoeDY8eOefWcomHSWnPs2DEcDkeNjpMagaib3lMgojnMnQJvjISr/wtN2lsdlWminMGcKiwir6AIh93mlXPGxcWRkpLCkSOBn0iF+RwOR8nIbE9JIhB112kUTJkPH0yEN/4AV34Mcb2tjsoULqcdMEYXeysR2O122rVr55VzCVEb8mhIeEdcElz3LYSEG4+JdvzP6ohM0Sk2kgm94wiS5/kigEgiEN4TczbcsAhiOsCHV8C6d6yOyOuS2jbmyQmJAbUmgRCSCIR3RTSDKV/BWUPgiz/DsscDbqyB1lrWLRYBRRKB8L7QSLjyI0i8EpY9Cl/eCevneG0xFSv9fiyLDg98w2fr91sdihBeI43Fwhw2O1zyIkS1gBVPw8/vgi5e69ePZy+NCA2moEiTbuXEc0J4mdQIhHmUMqancDY6kwRO89PZS6NKeg2ZOBW1ED4miUCYL+ek++1pKT4NwxvstiDCQmyk50qNQAQOSQTCfJXNUuqns5e6nHZr1yQQwsskEQjzuZu91BZS68VUrHblua0ZcHaM1WEI4TXSWCzMd7pB+PRiKja70aW0Wdeqj6un/jy8g9UhCOFVUiMQvpEwEe7abKxoNnUThMfAnCsh+7jVkdVYQWGRPBoSAUUSgfC9yOYw8V3IOAifXAeF/tUD577/bmL0syusDkMIr5FEIKwR3wcuehp2L4XFD1sdTY1EOYNlHIEIKNJGIKxzzp/g4AZY+Ry0SIQe462OyCMup52MvAIKizS2IJl8Tvg/qREIa13wGLTuZyx5eXCj1dF4JMphDCqTWoEIFJIIhLWCQ2Di2xDWGOZcBVnHrI6oWqfXJJBBZSJQSCIQ1otoBpPeNdY+/mRKvW887hHn4u6RHQkPlSerIjBIIhD1Q6veMPYZ+G05LKzfA806xkZyx/AOxETImgQiMJiWCJRSbyqlDiulNldS7lJKfamU2qCU2qKUutasWISf6HklnHczrHoBNsyxOppKFRZpDpzMkUdDImCYWSN4CxhVRfltwFatdSIwBHhaKRViYjzCH/zhEWgzwFjD4MDPVkfj1pGMPPrNWsL8DQetDkUIrzAtEWitlwNVDRvVQKRSSgERxfvW74fDwnw2u9F4HN4U5lwNmUesjqiCKKfRNiA1AhEorGwjeB7oAhwANgF3al1+0nrRIIXHwKT3IPsozJ0ChfXrC9dpt2G3KZlmQgQMKxPBBcB6oCXQE3heKRXlbkel1E1KqWSlVPKRI/XvL0RhgpY9Ydy/4ffv4X8PWB1NGUopohwyFbUIHFb2f7sWmKW11sBOpdRvQGdgdfkdtdavAq8CJCUlyarhDUXCRGPk8Y/PQ0Ee7FpszF7qijOmsLZwmUuX0y4DykTAsDIR7AWGAyuUUrFAJ2C3hfGI+mjEw/DrIlj31plt9WDN4zuGdyA6zG7JtYXwNtMSgVLqQ4zeQDFKqRTgIcAOoLV+GfgH8JZSahOggPu01kfNikf4KVswnEqvuP30mscWJYJLerWy5LpCmMG0RKC1vqKa8gPAH8y6vggg6ZV007RwzePD6bkcyzpFlxZum7WE8CsysljUf/VwzeNnF//KVa//ZNn1hfAmSQSi/nO35nGw09I1j6OKF7A3+joI4d8kEYj6L2EijH0OXPEYzUlA2wGW9xoqLNJknyq0LAYhvEUSgfAPpdc87n2t0ZV0/zrLwjk9FbWMJRCBQBKB8D8jZkB4M6MLqUVTVpcsTiPTTIgAIIlA+B9nNIx+Ag5tglUvWhJCr9bRPDOpJ82jHJZcXwhvkkQg/FOXcdBpNCx9FE7s8fnlW0Y7uaRXK6LDZMJc4f8kEQj/pBSMfhKCbDD/bvBx751TBUX8tPsYB07m+PS6QphBEoHwX6fnHNq1GDZ94tNLZ58qYNKrq1iw+ZBPryuEGSQRCP/W5wZjmcsF0yC7quUvvCvSIb2GROCQRCD8W5DNGGOQexK+/bvPLmsLUkSGBkuvIREQJBEI/9e8O/T7M6x/D3Z/57PLnh5dLIS/k0QgAsPg+6BRO5g/1ZiZ1AeinHbSc2R1VeH/JBGIwGB3wpjZcHw3LH/KJ5ecMbYrU0d08Mm1hDCTJAIRONoPhcQr4IdnIHWr6Zc776wmdG/lMv06QphNEoEILH/4J4RGGdNPFJk7Idz2Qxks3Jpq6jWE8AVJBCKwhDeBUY9ByhpIftPUS81N3sedc3429RpC+IIkAhF4EibBWUNh0cOQtt+0y7icdrJPFZJfWGTaNYTwBUkEIvAoZTQcFxXAN/eadpmo4qmo06ULqfBzkghEYGrcDoZMg1/mw7YvTbmErEkgAoUkAhG4zr8NYnvA1/dAbprXTx/lDAYkEQj/J4lABC6bHcY9CxmH4KM/wezuMCPa+Lnx4zqfvnfrxvz3lvPpGBtZ91iFsFCw1QEIYapWvaH9MGOG0tPS9hndS6FO6x67wuz0btO4jgEKYT2pEYjAd+SXitvyc2DxzDqdNq+gkLnJ+9h2ML1O5xHCapIIROBLP+B+e1pKnU6rNdzzyUaW/HK4TucRwmqSCETgc8XVbLuHHHYbIcFB0n1U+D1JBCLwDZ9uTEpXmt1pbK+jKIdd1iQQfk8SgQh8CRONxWtc8We2DZtep4bi01zOYOk+KvyeJALRMCRMhLs2w11bIdgJB9Z55bQuWZxGBADpPioaFlcrOP9WWPE09L0VWp1Tp9M9M6kXIcHy95Twb3IHi4an/50Q1gQWTje6/tRB6yZhNHc5vBSYENaQRCAaHocLBk+DPSvg12/rdKo1e47z+ordXgpMCGtIIhANU9K10Li9USsorP26w99tP8KjX29D17FmIYSVJBGIhslmhxEPGaOO179f69O4nHaKNGTmySL2wn9JIhANV5dxEH8eLH0UTmXV6hQyA6kIBJIIRMOlFIz8B2Qegh9fqNUpXCWL00iNQPgvSQSiYWt9HnQZCz88C5k1nzMoShanEQFAEoEQIx6GglxY9liNDz2ndSPWPDCCPm0bmRCYEL4hiUCIJu0h6TpY+zYc2VGjQx12G00jQwm2yf9Kwn/J3SsEwOD7wB4Gi2bU6LD8wiL+tXAHK3cdNScuIXzAtESglHpTKXVYKbW5in2GKKXWK6W2KKW+MysWIaoVHgMDpsL2r+D3lR4fZlOKfy/5lVW7jpkXmxAmM7NG8BYwqrJCpVQ08CIwTmvdDZhgYixCVK/vrRDZEr590OOpJ4KCFFEOmXhO+DfTEoHWejlwvIpdrgQ+1VrvLd5flnkS1goJg2EPwP61sGWex4dFyVTUws9Z2UbQEWiklFqmlFqrlPqThbEIYUi8App1g8UPQ0GeR4e4nHbSc2UcgfBfViaCYKA3cBFwAfB3pVRHdzsqpW5SSiUrpZKPHDniyxhFQxNkgz/MhBN7IPlNjw5xOe2yXKXwa1auR5ACHNVaZwFZSqnlQCJQof+e1vpV4FWApKQkmd1LmKv9cDhrCHz3uFFDcEZXufsbk/sQIt1HhR+z8u79HBiolApWSoUB5wHbLIxHCMPpqSdyTsL3s6vd3WG3ERSkzI9LCJOY2X30Q+BHoJNSKkUpdb1S6mal1M0AWuttwAJgI7AaeF1rXWlXUyF8qkUCJF4Oq16Ck/uq3HXJL6k8MG+TjwITwvvM7DV0hda6hdbarrWO01q/obV+WWv9cql9ntRad9Vad9daP2NWLELUytAHjJ9LHqlyt60H0nn/p73k5hf6ICghvE8ebApRmeh46HsLbPwIDm6odLeSGUhzpcFY+CdJBEJUZeDd4GxU5frGUSVTUUsiEP5JEoEQVXG4YPC9sHsZPHk2zIiG2d1h48clu5yZilrGEgj/ZGX3USH8g8MFKMgunlgubR98eYfxe8JEXE47TrtN2giE31L+tuh2UlKSTk5OtjoM0ZDM7m58+Zfnioe7NqO1RinpPirqN6XUWq11krsyeTQkRHXSUqrcLklA+DtJBEJUxxVX5faiIs1dH61n/sYDPgxKCO+RRCBEdYZPB7uz7LZgp7EdYyrqbzYfZMO+k76PTQgvkEQgRHUSJsLY54w2AYofA3W80NhezJh4TnoNCf8kvYaE8ETCxDNf/B9Mgt2LjbmIiiekczllcRrhv6RGIERNDXsQctNg5b9LNskqZcKfSSIQoqaa94DufzQmpMs0FtaLa+QkPNRmcWBC1I4kAiFqY+gDUJALK54G4JnLe/H65D4WByVE7UgiEKI2mrSHXlcbq5id3Gt1NELUiSQCIWpr8L2Agu8e58sNB7jq9VUUFvnXSH0hQBKBELXnioM+N8D6D8g7tJ0fdh4jUxaxF36o0kSglJKupUJUZ+DdYA/jvD0vAbImgfBPVdUIVvssCiH8VXgM9L2V+IP/o5vaI11IhV+qKhHITFpCeKLf7eSHRPPX4I8kEQi/VNXjn6ZKqbsrK9Ra/8uEeITwPw4Xab1vY+iP/2Tz4TVw9oVWRyREjVRVI7ABEUBkJS8hRLGYobdDRHO6//JspUtaClFfVVUjOKi1numzSITwZyFhMPge+OovsHMxdBhhdURCeEzaCITwkivWdiTD2QoWPwxFRVaHI4THqkoEw30WhRAB4JfDufyv6bVwaCNs+9zqcITwWKWJQGt93JeBCOHvopx2ljuGQtMusOSfUCiDy4R/kJHFQniJy2knLbfImKb62K+wcY7VIQnhEUkEQnhJyeI0nS+ClufAsllQkGd1WEJUSxKBEF7SvZWLDs0iQCljPeO0fbD2LavDEqJaMp+QEF5y36jOZ96cNQTaDoTlTxrTVYeEWxaXENWRGoEQZlAKhj8EWUeMlcyEqMckEQjhJR8n72PgE0vIzS80NsT3gU6j4YfnIOeEtcEJUQVJBEJ4SX5hEfuO55SdeG7oA5CXbiQDIeopaSMQwkuiHHYA0nLyiY1yGBubd4ce42Hlv2HDHMg4aCxoM3w6JEy0MFohzpAagRBe4nIaiSC9/FTULRKhKB8yDgDa6E305R2w8WPfBymEG5IIhPCS04mgwpoEP71Scef8HFgsczqK+kESgRBe0jQylBFdmpUkhBJpKe4PqGy7ED4mbQRCeEnLaCevT+5TscAVZzwOcrddiHpAagRCmG34dLA7y26zO43tQtQDkgiE8KKhTy1j1je/lN2YMBHGPgeu+DPb+twovYZEvSGJQAgvyssv5Gimm4nmEibCXZvhwSPQ+CzYuQiKCn0foBBuSCIQwouiTs9AWpngEBgxAw5vhfXv+ywuIapiWiJQSr2plDqslNpczX59lFKFSqnxZsUihK9EOe0VxxGU12UcxJ0LSx+FU1m+CUyIKphZI3gLGFXVDkopG/A48D8T4xDCZ1zV1QjAmJDuD48Yo4x/fNE3gQlRBdMSgdZ6OVDdcpd/Bv4LHDYrDiF8aVCHGIZ1blb9jq3Pgy5j4YdnIFNuf2Ety9oIlFKtgEuBl62KQQhvu+b8ttxbel2CqgyfAQW58N3jpsYkRHWsbCx+BrhPa11t1wml1E1KqWSlVPKRI0fMj0yIOigq0mitq98x5mzofS0k/weO/mp+YEJUwspEkATMUUrtAcYDLyqlLnG3o9b6Va11ktY6qWnTpj4MUYiaeXfV77R/4GtOZFfTTnDa4PvAHgaLZpgalxBVsSwRaK3baa3baq3bAp8At2qtP7MqHiG8IcxuQ2s3M5BWJqIpDLgTfpkPv/9obnBCVMLM7qMfAj8CnZRSKUqp65VSNyulbjbrmkJYrdIZSKvS9zaIbAEL/w6ePFISwstMm3ROa31FDfadYlYcQvhS1Ok1CXJrkAhCwoyVzL64HbZ+Bt0uNSc4ISohI4uF8KJa1QgAel4JzbrBooeh4JQJkQlROUkEQnhRs8hQpvRrS+vGYTU7MMgGI2fCid8g+U1zghOiEpIIhPCiRuEhzBjXjYS46JoffPZwaDfYGFeQc9LboQlRKUkEQnhZbn4hWXkFNT9QKfjDPyDnhDHiWAgfkUQghJcNeHwJ//x6W+0ObpEICZNg1UuylKXwGUkEQnhZlMODGUirMuwBoxvpkke8F5QQVZBEIISXVbsmQXWiW0Pfm2HDHDi40XuBCVEJSQRCeJlHaxJUZ8Dd4IyGhbKusTCfJAIhvMzltJOeW4vG4tKc0TDoXti91FjWUggTSSIQwsvGJLRg8vlt6n6iPjdAo7aw8CFZ31iYShKBEF52QbfmTOnfru4nCg6B4Q9B6majvUAIk0giEMLLcvML2Xc8m4LCorqfrNul0Ko3LLgf/tUNZkTD7O6w8eO6n1uIYpIIhPCyz37ez8AnlnI4I6/uJ1MK2g+DvJOQngJoSNsHX94hyUB4jSQCIbysVjOQVsXdY6H8HFg80zvnFw2eJAIhvKxkBlJPVymrTmUjjGXksfASSQRCeFmUo5ZTUVfGFVez7ULUkCQCIbzMVfJoqI5jCU4bPh3szrLbgp3GdiG8QBKBEF7WNDKU6WO6khjn8s4JEybC2OfAFQ8oY1v7ocZ2IbzAtKUqhWionCE2rhvghXEEpSVMPPPFP+8W2PQxpG6B2G7evY5okKRGIIQJdh3JZN/xbHNOfsE/weGCL+6QEcfCKyQRCGGCP72xmtkLd5hz8rDGMGoW7E+WZS2FV0giEMIEdZ6Kujo9JhgDzRY9DGn7zbuOaBAkEQhhApcz2HsDytxRCi76FxQVwDf3mncd0SBIIhDCBC6zawQAjdvB0L/BL/Nh25fmXksENEkEQpjAWK7SS+MIqtL3NmjeA76+B3LTzL+eCEiSCIQwweXntmbGOB907bQFw9hnITNV5h4StSaJQAgT9G7TiFHdm/vmYq16w3k3w5o3YO9PvrmmCCiSCIQwwdHMPH7YeZTcfB/18x/6gDH30Jd3QsEp31xTBAxJBEKY4Ptfj3LV6z9x4GSOby4YGmH0IjqyDVY+65trioAhiUAIE0Q5jdlbTO85VFrHP0C3y+C7J+HoTt9dV/g9SQRCmMDrM5B6atQssDtg/lTQ2rfXFn5LEoEQJihZnMaXNQKAyFgY+Q/YswJ+fs+31xZ+SxKBECY4vThNuq8TAUCva6BNf/j2Qcg87PvrC78jiUAIEzQKD+G1PyUxtHMz3188KAjGPAP52bDgb76/vvA7kgiEMIHdFsTIrrG0inZWv7MZmnaEgX+FzZ/ArwutiUH4DUkEQpjkh51HWbf3hHUBDJgKMZ1g/t1wKsu6OES9J4lACJM89MUWXl+x27oAgkON6SfS9sLSR62LQ9R7slSlECbxyQyk1WlzPvS+Fn58HjbNNRqPXXHGwvey5rEoJjUCIUwS5Qj2zQyk1WnR0/iZmQpoSNsHX94BGz+2MipRj0giEMIk9aJGALDiqYrb8nNktlJRwrREoJR6Uyl1WCm1uZLyq5RSG4tfK5VSiWbFIoQVTF+u0lNpKTXbLhocM2sEbwGjqij/DRistU4A/gG8amIsQvjc9QPa8e7151odhtEmUJPtosExLRForZcDx6soX6m1Pt23bhUgd6UIKG2ahJMQF211GEbDsL38eAZljDMQgvrTRnA98I3VQQjhTfuOZzNn9V7rHw8lTISxz4ErHlAQXjzaefvXUFRkaWiifrC8+6hSaihGIhhQxT43ATcBtG7d2keRCVE3Ww6kMe3TTfSIc+FyuqwNJmFi2e6ia16Hr/4Cyx6FYQ9aF5eoFyytESilEoDXgYu11scq209r/arWOklrndS0aVPfBShEHURZNQOpJ5KuNyanW/4kbP3c6miExSxLBEqp1sCnwDVa6x1WxSGEWc7MQFoPxhKUpxRc9DTE9YF5t0DqVqsjEhYys/voh8CPQCelVIpS6nql1M1KqZuLd5kONAFeVEqtV0olmxWLEFYoWZymPtYIwJiCYuK7xjKXc66EHAvnRRKWMq2NQGt9RTXlNwA3mHV9IazmCqvHj4ZOi2phJIO3LoJProer5kKQzeqohI/Vl15DQgSciJBgvr1rEBP7xFsdStVanwcXPQW7FsOSf1gdjbCA5b2GhAhUQUGKjrGRVofhmd5T4MB6+H42NE+A7pdZHZHwIakRCGGiT9elsGDzIavD8MyFT0B8X/j8NjjkdmYYEaAkEQhhojd/+I2Pk/dZHYZngkNg4jvgcBmNx9mVTgwgAowkAiFMVG9mIPVUZCxMeg8yDsIn10FhPez6KrxOEoEQJopy+FkiAIhLgov+BbuXwuIZVkcjfEAai4Uwkctpr7/jCKpyzjVwcAOs/Dc0T4SECVZHJEwkNQIhTOR3j4ZKG/UYtOkPX/zZSAoiYCmttdUx1EhSUpJOTpZByMI/pOXkU1SkaRQeYnUotZN5BF4dbKxoZndA+kFZ89hPKaXWaq2T3JVJjUAIE7mcdv9NAgARTY0xBjnHIf0AsuZxYJJEIISJdh7O5IkFv3A4PdfqUGpv3TsVt8maxwFFEoEQJko5kc2Ly3ax70SO1aHUnqx5HPAkEQhhono/A6knKlvbOLK5b+MQppFEIISJ6vXiNJ5yu+YxcCoHUrf4Ph7hdZIIhDBRSY0g148TQfk1j13xMPwhCHHCm6PgtxVWRyjqSAaUCWGiM6uU+XEigIprHgP0mADvj4f3LoNLX4buf7QmNlFnkgiEMFFIcBDbZo7CYQ/Aynd0PFy3AD680piXKOMQnH+b1VGJWgjAu1OI+sUZYkMpZXUY5nA2gmvmQZdx8L/74X8PQFGR1VGJGpJEIITJXv5uF2+v3GN1GOaxO2DCW3Du/8GPz8OnN0BBntVRiRqQRCCEyRZtTeWbzQetDsNcQTa48HEYORM2/xfe+yPkplkdlfCQJAIhTGbMQNoA5vVXCvrfCZe+Cnt/hDcvLJ6WQtR3kgiEMJlfz0BaG4mT4Kq5cPJ3eH0kHN5mdUSiGpIIhDBZlNPu3+MIaqP9MLj2ayjKhzcvgCWPwOzuMCPa+CkT1tUrkgiEMJnLaUcBRUX+NeV7nbVIhOsXgs0By580Zi2V2UvrJUkEQphs6ogObJxxAUFBAdqFtCqN2oDNzXAlmb20XpFEIITJAnYMgacqazCW2UvrDUkEQphs64F0bv9gHb8dzbI6FGtUNntpUDD8vtK3sQi3JBEIYbL03HzmbzzIgZN+vCZBXbibvdQWAiHh8J8LYe61cHKfNbEJQBKBEKYLiDUJ6sLd7KUXvwB3b4PB02D7N/B8Eix9FE410FqTxWTSOSFM5gqENQnqyt3spQBD/wa9roZFD8F3j8O6d43RyT3GGwPUhE9IjUAIkwXE4jRmio6H8W/CtQsgopkxV9GbF8D+tVZH1mBIIhDCZOEhNmKjQrE1xO6jNdHmfLhxqfHY6Phv8Now+OxWY3rrjR/LgDQTKa39a5BLUlKSTk5OtjoMIYSZctNhxdOw6kVAgS4yRimfZnca7Q7uHjcJt5RSa7XWSe7KpEYghKh/HFEw8mG4dRWooLJJAGRAmpdJIhDCBx77ehuPfS2Tr9VYk/ZQkOu+LG2f8dhI1Jn0GhLCB7YeTCcjtwFMRW0GV1zxPEVuPN0Z2vSHbpdA14uNxmZRY1IjEMIHGuQMpN7ibkCa3QkjHoYh0yDrCHz9V3i6E7w9FpLfhKyjZfeXxuYqSY1ACB8wFqeRRFArpxuEF8805idyxRnJ4fT2wfcZax5smQdbPoX5d8FXf4V2A6HbZaALjfWU84tHdp+e/bT0uRs4SQRC+ECUw1icRmstk9DVRmUD0sAYeBbb1XgNvR9St5xJCqe/8Ms73dgsiQCQR0NC+ER8YycdmkVyqrDI6lACm1LQvDsM/zv8eR383/LK903bB6tegt+WQ9axqs8b4I+WZByBECKwze7uvrFZBRnjE06LiIXYbtCs65mfTTvDti+MmkV+qUkD/XAcQ1XjCEx7NKSUehMYAxzWWnd3U66AZ4HRQDYwRWu9zqx4hBAN1PDplX+RtxsMh7dA6lbjkdLhLbD6NSjMM/ZTNqOWUVSux1d+DiycbvRUCg6tPoaNH1fexuGJuh5fDTPbCN4CngfeqaT8QqBD8es84KXin0IEnM3703hg3iYeuaQHPeJcVofTsFTX2BwZa6yxfFphARzfXZwgthjLbLqTcRAeaQYOF4Q3LX7FlPq9+P2hzfDj82fGQ9S0sXrjx2UTmQmN3aYlAq31cqVU2yp2uRh4RxvPplYppaKVUi201gfNikkIqxQWaTakpLFw6yEOpecSbFMM7WT0ed+w7ySHM/LK7O+wBzGwQ1MA1v5+guNZp8qUh4fa6Nc+BoDVvx2vMKGdy2nn3HaNAVi56yhZeYVlyhuHh9C7TSMAvv/1KDn5ZcubRYaSGB8NwLLth8kvLPsIuYXLQfdWRkJbvC2V8ssxxzVy0qVFFEVFmsW/HK7webRtEkaH2EhOFRTx3Y4jFcrbNw3nrKYR5Jwq5PudRyuUd4qNpHWTMDJy81m1+3iF8q4to2gV7eRk9inW7DkBwYPhgsUl5QlxLmKBo5l5/Lz3ZIXje7VuQ0y3jqTGX0hU8vs4s92ssuZsRHrPmzh55AAheccIyT1OyMlfCMn7AXveCRRVPHbPz4F5N3Nq2VNk46Aw2ElhcDiFwWEUBofRvGkTgh2RHD1lIyr5BULycyoe78XGbit7DbUCSj+4SyneViERKKVuAm4CaN26tU+CE8KbmrscKAXPLdkJQHSYnfXT/wDAy9/t4pvNZUfItop28sM046/UZxbtYMWvZb8MO8ZG8O1dgwGY9c021pX7MuvVOpp5t/YH4OEvtrI9NaNM+cAOMbx7vVEBn/bpRlJOlP2iGdWtOS9f0xuAuz5az4nssonmsnNa8a+JPQG45b11FRrB/3R+G2Ze3J2CIs2N71Rs07tlSHvuG9WZrLwCt+X3XNCJ24aezbGsPLflM8Z2ZUr/dhw4meu2/InxCUxMimfXkSy35S9edQ6je7Rg64F0t+VvX3cugzs2Zd3vJ/jm5CXMsr9OmDqTjAttTmwXPsG3+f3469INFY5fcEc/Okfl8/kPGxi7cjxu5xvUhewNimN/6hHC1HHCOUgYuYSrXGz7CiA/ixg3h5Xw4lKfpjYWF9cI5lfSRvAV8JjW+vvi94uBe7XWVc49K43Fwl/tO55d8pd7kFJ0bRkFwN5j2RUGm9ltQXRqHgnAb0ezyMor+4w6NDiIDrFG+a4jmeScKvsXvTPERvumEQD8mppBXkHZL+rw0GDaxYQDsP1QBvnlvsijHHZaNwkDjKU2i8p9T7icduIbG+Wb96dV+Lc2Dg+hZbSToiLN1oPpFcpjIkJp7nJQUFjEL4cyKpQ3iwqlWaSDvIJCfk3NrFDe3OUgJiKU3PxCdh6uWN4q2kmj8BCy8grcLhEa3ygMV5idjNx8fj+WXaG8TZMwIh120rLz2XciG9fOecSueQJ75gHyI1pSNGw6jnMu50TWKfa7WXnu7GYROOw2jmbmEfVyL0Iy91fYB1c8h29I5nB6XoWizs0jCVaQevwEjd7sR0i2m6k0XPFw1+aK2ytRVWOxlYngFWCZ1vrD4vfbgSHVPRqSRCCE8Cvln/FDzXod1fX4YvV19tEvgD8pQ18gTdoHhBABx91SnTX5Eq/r8R4wrUaglPoQGALEAKnAQ4AdQGv9cnH30eeBURjdR6/VWlf7p77UCIQQouYsGUegtb6imnIN3GbW9YUQQnhGppgQQogGThKBEEI0cJIIhBCigZNEIIQQDZwkAiGEaOAkEQghRAMniUAIIRo4v1uYRil1BPi9lofHABWnMrRefY0L6m9sElfNSFw1E4hxtdFaN3VX4HeJoC6UUsmVjayzUn2NC+pvbBJXzUhcNdPQ4pJHQ0II0cBJIhBCiAauoSWCV60OoBL1NS6ov7FJXDUjcdVMg4qrQbURCCGEqKih1QiEEEKUEzCJQCk1Sim1XSm1Uyk1zU25Uko9V1y+USl1jqfHmhzXVcXxbFRKrVRKJZYq26OU2qSUWq+U8uoiDB7ENUQplVZ87fVKqemeHmtyXPeUimmzUqpQKdW4uMzMz+tNpdRhpZTbtQEtvL+qi8uq+6u6uKy6v6qLy+f3l1IqXim1VCm1TSm1RSl1p5t9zL2/tNZ+/wJswC7gLCAE2AB0LbfPaOAbQAF9gZ88PdbkuPoBjYp/v/B0XMXv9wAxFn1eQzCWGa3xsWbGVW7/scASsz+v4nMPAs4BNldS7vP7y8O4fH5/eRiXz+8vT+Ky4v4CWgDnFP8eCezw9fdXoNQIzgV2aq13a61PAXOAi8vtczHwjjasAqKVUi08PNa0uLTWK7XWJ4rfrgLivHTtOsVl0rHePvcVwIdeunaVtNbLgeNV7GLF/VVtXBbdX558XpWx9PMqxyf3l9b6oNZ6XfHvGcA2oFW53Uy9vwIlEbQC9pV6n0LFD7KyfTw51sy4SrseI+ufpoFvlVJrlVI3eSmmmsR1vlJqg1LqG6VUtxoea2ZcKKXCMJY5/W+pzWZ9Xp6w4v6qKV/dX57y9f3lMavuL6VUW6AX8FO5IlPvL9OWqvQx5WZb+e5Qle3jybG15fG5lVJDMf5HHVBqc3+t9QGlVDNgoVLql+K/aHwR1zqMIemZSqnRwGdABw+PNTOu08YCP2itS/91Z9bn5Qkr7i+P+fj+8oQV91dN+Pz+UkpFYCSeqVrr9PLFbg7x2v0VKDWCFCC+1Ps44ICH+3hyrJlxoZRKAF4HLtZaHzu9XWt9oPjnYWAeRjXQJ3FprdO11pnFv38N2JVSMZ4ca2ZcpVxOuWq7iZ+XJ6y4vzxiwf1VLYvur5rw6f2llLJjJIH3tdafutnF3PvL2w0fVrwwaja7gXacaTDpVm6fiyjb2LLa02NNjqs1sBPoV257OBBZ6veVwCgfxtWcM+NMzgX2Fn92ln5exfu5MJ7zhvvi8yp1jbZU3vjp8/vLw7h8fn95GJfP7y9P4rLi/ir+d78DPFPFPqbeXwHxaEhrXaCUuh34H0Yr+pta6y1KqZuLy18GvsZoed8JZAPXVnWsD+OaDjQBXlRKARRoY1KpWGBe8bZg4AOt9QIfxjUeuEUpVQDkAJdr486z+vMCuBT4VmudVepw0z4vAKXUhxg9XWKUUinAQ4C9VFw+v788jMvn95eHcfn8/vIwLvD9/dUfuAbYpJRaX7ztfowk7pP7S0YWCyFEAxcobQRCCCFqSRKBEEI0cJIIhBCigZNEIIQQDZwkAiGEaOAkEQghRAMniUAIIRo4SQRC1JFSqk/xHPEOpVR48Zzy3a2OSwhPyYAyIbxAKfUI4ACcQIrW+jGLQxLCY5IIhPACpVQIsAbIxZjXp9DikITwmDwaEsI7GgMRGCtMOSyORYgakRqBEF6glPoCY3WodkALrfXtFockhMcCYvZRIayklPoTxqyeHyilbMBKpdQwrfUSq2MTwhNSIxBCiAZO2giEEKKBk0QghBANnCQCIYRo4CQRCCFEAyeJQAghGjhJBEII0cBJIhBCiAZOEoEQQjRw/w8ciHRTkQ9+igAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#Primero definimos los parámetros del problema\n", "\n", "L = 2.0 #longitud del dominio en x\n", "N = 20 #numero de nodos espaciales en x\n", "dx = L/(N-1)\n", "tf = 0.2 #tiempo final que queremos simular\n", "b = 0.3 # valor de la conductividad térmica\n", "sigma = 0.3 #sigma es un parametro que controla la estabilidad del metodo numérico\n", "\n", "#después calculamos el paso de tiempo y el numero de pasos necesario\n", "dt = sigma*dx**2/b #se calcula el paso de tiempo para la integración temporal\n", "\n", "nt=int(tf/dt) # número de pasos de tiempo necesarios\n", "\n", "\n", "print(\"La longitud del dominio es:\",L)\n", "print(\"El numero de nodos espaciales es:\",N)\n", "print(\"El numero de pasos temporales es:\",nt)\n", "print(\"El paso de tiempo es es:\",dt)\n", "print(\"El tiempo final es:\",(nt-1)*dt)\n", "print()\n", "\n", "#inicializamos los arrays\n", "x = np.zeros(N)\n", "t = np.zeros(nt)\n", "T = np.zeros((nt, N))\n", "\n", "#damos valores a las posiciones espaciales\n", "for i in range(N):\n", " x[i]=dx*i\n", "\n", "#damos valores a los tiempos\n", "for i in range(nt):\n", " t[i]=dt*i\n", "\n", "#definimos la condición inicial\n", "def initial_condition1(x):\n", " return np.piecewise(x, [x < 1, x >= 1], [2, 1])\n", "\n", "#imponemos la condición inicial en el tiempo \"cero\"\n", "T[0,:] = initial_condition1(x)\n", "\n", "#calculo de evolución de la temperatura en el tiempo\n", "for n in range(0,nt-1): #bucle temporal \n", " for i in range(1, N-1): #bucle espacial entre 1 y N-2 (sin incluir primer y ultimo nodo)\n", " T[n + 1, i] = T[n,i] + b*dt/dx**2 * (T[n,i+1] + T[n,i-1] - 2*T[n,i])\n", " #condiciones de contorno en primer y ultimo nodo \n", " T[n + 1, 0]=2\n", " T[n + 1, N-1]=1\n", "\n", "n=nt-1 #instante de tiempo en el que pinto la solución\n", "\n", "fig, ax1 = plt.subplots(figsize=(6, 5))\n", "ax1.plot(x,T[0, :],'--',label='Condicion inicial')\n", "ax1.plot(x,T[n,:],'-o',label=f'Solución en t= {(n)*dt:.3f} s')\n", "ax1.set_xlabel(\"x\") \n", "ax1.set_ylabel(\"T\") \n", "ax1.legend()" ] }, { "cell_type": "markdown", "metadata": { "id": "W5lkkhCqVTOu" }, "source": [ "Ahora vamos a hacer una representación de la temperatura en función del espacio y del tiempo." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "e51Qwe_dVTDN", "outputId": "beb4f3c2-0d6e-4822-eeee-31bd9776e7fb" }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAj30lEQVR4nO3df9RdVX3n8fen+VEkhELMENIQBDUq0fJrMmCLo0GqJHQ0pFMrqQtZDE7EmlFmaldZrlni1Fmr1EotTpH4oBlgLX4MHUlNayQwqKUW0QQICSGgMfLjISExRA1WIYZ854+zHzy53B/nPs8995577+e11l333rP3PmffJyf3e8/e++ytiMDMzOzXel0BMzOrBgcEMzMDHBDMzCxxQDAzM8ABwczMEgcEMzMDHBBswEmaL2lDh/Z1sqR7O7EvsypyQLC+JelnucdBSb/IvX9fyvYp4DPj3H9Ieu3Y+4jYBPxE0rs6UH2zynFAsL4VEUeMPYAngXfltt0kaTZwNvD3HTzsTcAHO7g/s8pwQLBB9g7ggYh4HkDSayTtlXR6ev+bkvZIWlhbUNI96eVD6Yrjven9N4FzJP166bU36zIHBBtkvwU8NvYmIn4A/Blwk6TDgf8NXB8R36wtGBFvTS9PSVcc/ydtfxr4JfD6kutu1nWTe10BsxIdBTyb3xAR16U+gO8AAbx7HPt9Lu3bbKD4CsEG2Y+B6XW2Xwe8CfhfEfHCOPY7HfjJBOplVkkOCDbINgGvy2+QdATwN8CXgE9KmtHODiX9JjCVXFOU2aBwQLBBdhdwuqTDctuuBu6PiA8AXwVWjiVI+qSkb+by7gJeXbPPhcDXx3llYVZpDgg2sCJiF/B1YAmApCXAIuDSlOW/kQWMsXsW5gL/ktvFJ4EbJP1E0h+mbe8jF0TMBom8QI4NMknzgRuAM6LFyS5pI3BORDzbIP23gJGI+O2OV9SsAhwQzMwMcJORmVllSZor6RuStkraIumjdfK8QdK3Jb0g6WM1aYskPSZpm6TLWx7PVwhmZtWUpl+ZHREPSJoO3A+cHxGP5PIcA7wKOB/4cUR8Jm2fBHyP7I79UWA9sCxftpavEMzMKioidkbEA+n1c8BWYE5Nnt0RsZ7sDvq8M4BtEbE9IvYDt5IGWDQyFHcqHz3j12LOcUPxUQH48YuH97oKh9h34LDWmdq0f//E/z21Xx2oCUwqYQDqpBeqdeWu5/f3ugpds+/Aj/ZExL+ZyD7esvCw+Mnegy3zbdn8yy3A87lNIxExUi+vpBOA08jusi9iDvBU7v0ocGazAkPxLTnnuMnc9tUJ/fv2ldX7Tut1FQ6x7pmTOr7PJ5+aOeF9/PqTUztQE5j+ROe/vH9je7Vuc5j66Givq9A1dzzz+Scmuo+f7D1Y6DvnjcfveD4iFrTKl26o/DJwWUTsK1iNer94mp6sbjIyM6swSVPIgsFNEXF7G0VHye6tGXMcsKNZAQcEs4rx1YGNkSSyaVa2RsRft1l8PTBP0omSpgIXAGuaFRiKJqNhs/TIByvXbGRm43IWcCGwOd04CfBx4HiAiFgp6VhgA3AkcFDSZcD8iNgnaQWwDpgErIqILc0O5oBgpTv32K2l9COYDbqI+Bb1+wLyeZ4haw6ql7YWWFv0eG4yMjMzwAHBzMwSBwQzMwMcEMysCY8wGi4OCANq6ZEP9roKZtZnHBDMzAxwQDCbkDKmrTDrlVIDQqu5uBvN4y3p9ZI25h770s0WY+vePp1LO6/Mz2DWTVW7S9mGS2k3pqW5uK8hNxe3pDU1c3HvBT5CNo/3SyLiMeDU3H6eBlbnsnx2bM5vMzPrjDKvEFrOxd1kHu+8c4AfRMSEZyC03jn32K29roKZtVBmQKg3F/ecBnmbuQC4pWbbCkmbJK2SdHS9QpKWS9ogacPeAvOSDyKPNLKJ8JDT4VNmQGh7Lu6X7SCboe/dwN/lNl8LvIasSWkncFW9shExEhELImLBjBnuOzcza6XMb8q25+KuYzHwQETsGtsQEbsi4sWIOAhcR9Y0ZWZmE1RmQGh7Lu46llHTXJQWnR6zFHh4QrW0vnT83D29roLZwCltlFFEHKg3F7ekS1N6q3m8DycbofTBml1/WtKpZM1Pj9dJNzOzcSh1PYR6c3FHxMrc62bzeP8ceGWd7Rd2uJpmVsMdysPJva0DziONzKwoBwSzivBdytZrDghmZhWW7rfaLanuABpJR0tane7N+q6kN+XSHpe0OU3zs6HVsRwQrGt8t7LZuFwPLGqS/nFgY0ScDLwfuLom/eyIODUiFrQ6kAOCmVmFRcQ9ZPO+NTIfuDvlfRQ4QdKs8RzLAcHMDuERRn3nIeD3ASSdAbyKX43eDOBOSfdLWt5qR6UOO7VqWHrkg6zed1qvq2E2NH784uEF/8/tmFnTtj8SESNtHu5K4GpJG4HNwIPAgZR2VkTskHQMcJekR9MVR10OCGZmvbOnSNt+MxGxD7gYQJKAH6YHEbEjPe+WtJpsqp+GAcFNRmZmfUzSUWl6IIAPAPek2R6mSZqe8kwD3kmLqX58hWBmVmGSbgEWAjMljQJXAFPgpZkfTgJulPQi8AhwSSo6C1idXTQwGbg5Iu5odiwHBDN7iTuUqycilrVI/zYwr8727cAp7RzLTUZDwlNYVJvvUrYqcEAwMzPAAcG6zHcrm1WXA4KZmQEOCGbjNv2JtpYIrzx3KJsDwhBxx7KZNeOAYGZmgAOCmZklpQYESYskPSZpm6TL66S/QdK3Jb0g6WM1aXUXdpA0Q9Jdkr6fno8u8zOYmQ2L0gKCpEnANcBisvm6l0maX5NtL/AR4DMNdlNvYYfLgbsjYh7ZHOAvCzTWmPsRrB53KBuUe4VwBrAtIrZHxH7gVmBJPkNE7I6I9cAv29jvEuCG9PoG4PwO1NX60PFz9/S6CmYDpcyAMAd4Kvd+NG0rqtHCDrMiYidAej5mwjU1M7NSJ7dTnW3tDNxua2GHlx08CyLLAWbPmdTGYc26y/MYWVWUeYUwCszNvT8O2FG0cH5hB2BsYQeAXZJmA6Tn3Q3Kj0TEgohYMGOGB1OZNeL+AxtT5jflemCepBPT4g0XAGuKFGyxsMMa4KL0+iLgKx2t9RDodcey5zMyq6bSmowi4oCkFcA6YBKwKiK2SLo0pa+UdCywATgSOCjpMrIRSTNpvLDDlcBtki4BngTeU9ZnMDMbJqUukBMRa4G1NdtW5l4/Q9aUVGsfDRZ2iIhngXM6WE0zM8N3KpuZWeKAMKR63Y9g1eAOZctzQDAzqzBJqyTtlvRwg/TfkPQPkh6StEXSxbm0ptMH1XJAMDOrtuuBRU3SPww8EhGnAAuBqyRNLTh90CEcEMzMKizdkLu3WRZgurJhmUekvAcoMH1QrVJHGVm1LT3yQVbvO63X1bAecf9BefYdOIx1z5xUIOdXZ+ZncwZGImKkzcP9Ldn9WTuA6cB7I+KgpHrTB53ZbEcOCGY95Gkrht6emtmcx+NcYCPwduA1ZFP9/DPjmD7ITUZmZv3tYuD2yGwDfgi8gXFMH+SAYD3h6SvMOuZJ0s26kmYBrwe2M47pg9xkZDaE3H/QPyTdQjZ6aKakUeAKYAq8NPPDp4DrJW0mayb6s4jYk8q+bPqgZsdyQBhy7lgen+lPtDOTu9n4RcSyFuk7yCYArZf2sumDmnGTkZmZAQ4IZmaWOCCYmRnggGBmZokDgvW14+fu6XUVzAaGA4LZkPGQU2vEAcG8NkKPeNoKqxoHBDMzAxwQrIc8fYVZtZQaEFqt1iPpDZK+LekFSR/LbZ8r6RuStqYVgD6aS/ukpKclbUyP88r8DGZmw6K0qStyq/W8g2zWvfWS1kTEI7lse4GPAOfXFD8A/ElEPCBpOnC/pLtyZT8bEZ8pq+5mZsOozCuElqv1RMTuiFgP/LJm+86IeCC9fg7YCswpsa5mQ8EjjKyZMgNCvdV62v5Sl3QCcBrwndzmFZI2pcWnj55QLQ3wSCMzKzcgtL1az8t2IB0BfBm4LCL2pc3Xkq0KdCqwE7iqQdnlkjZI2rB378F2DmtmNpTKDAhtr9aTJ2kKWTC4KSJuH9seEbsi4sWIOAhcR9Y09TIRMRIRCyJiwYwZHkxlZtZKmd+Uba/WM0aSgC8BWyPir2vSZufeLgUe7lB9zcyGWmmjjCLiQL3VeiRdmtJXSjoW2AAcCRyUdBkwHzgZuBDYLGlj2uXH02IPn5Z0Klnz0+PAB8v6DGZl6cVdyu5QtlZKXTGt3mo9acm3sdfPkDUl1foW9fsgiIgLO1lH+xWvnmY23Ny4bj3lu5XNmkujKXdLatg8LmlhulF3i6R/ym1/XNLmlLah1bG8prKZWbVdD/wtcGO9RElHAZ8HFkXEk5KOqclydkQUmifeVwhmbZr+RFujp80mJCLuIZvVoZE/Am6PiCdT/t3jPZavEMzMOmz//sk8+dTMIlln1jTljETESJuHex0wRdI3genA1RExdjURwJ2SAvhCq307INgh3LE8mDzCqLL2RMSCCe5jMvBvgXOAVwDflnRfRHwPOCsidqRmpLskPZquOOpyk5GZWX8bBe6IiH9NfQX3AKcARMSO9LwbWE2DG3nHOCCYmfW3rwD/XtJkSYcDZwJbJU1Ls0UjaRrwTlrcyOsmI+t7x8/dU7S91qzvSLoFWEjW3zAKXAFMgey+rojYKukOYBNwEPhiRDws6dXA6mziByYDN0fEHc2O5YBg1mVeS9naERHLCuT5K+CvarZtJzUdFeUmI7MB5w5lK8oBwV6m22sj+G5ls2pwQDAzM8ABwczMEgcEMzMDHBCsAa+xPBjcoWztKDTsNK149rr09rGI+GV5VTIzs15oGRAkLQRuIFudTMBcSRc1mw/DzMz6T5ErhKuAd0bEYwCSXgfcQjaZkpm1wTelWZUV6UOYMhYMANIMelPKq5JVhfsR+pv7D6xdRQLCBklfSku0LZR0HXB/2RWz4eKb08x6r0hA+BCwBfgI8FHgEeDSIjuXtEjSY5K2Sbq8TvobJH1b0guSPlakrKQZku6S9P30fHSRupiZWXMtA0JEvEC2nuf/AD4BXJO2NSVpEnANsBiYDyyTNL8m216yQPOZNspeDtwdEfOAu9N7MzOboJYBQdLvAT8AriYLDNskLS6w7zOAbRGxPSL2A7cCS/IZImJ3RKwHaoexNiu7hGzUE+n5/AJ1MeuIfllP2f0HNh5FmoyuAs6OiIUR8TbgbOCzBcrNAZ7KvR9N24poVnZWROwESM/H1NuBpOWSNkjasHfvwYKHtVruWDYbHkUCwu6I2JZ7vx3YXaCc6mwr+vNqImWzzBEjEbEgIhbMmOEbss3MWilyH8IWSWuB28i+lN8DrJf0+wARcXuDcqPA3Nz744AdBevVrOwuSbMjYqek2RQLTmZm1kKRn86HAbuAt5Et4/YjYAbwLuA/NCm3Hpgn6cQ09cUFwJqC9WpWdg1wUXp9Edl6olYiNxt1RrduSnP/gY1XyyuEiLh4PDuOiAOSVgDrgEnAqojYIunSlL5S0rHABuBI4KCky4D5EbGvXtm06yuB2yRdAjxJdsViQ87rKtugkrSK7Mf37oh4U530hWQ/jH+YNt0eEX+e0haRDQiaRLbW8pXNjlVkLqMTgf8CnJDPHxHvblU2ItYCa2u2rcy9foasOahQ2bT9WeCcVse2/nPusVtZ98xJva6GWdVcTzbC88Ymef45Ig5psckN338HWTP8eklrIuKRRjsp0ofw98CXgH8APFzHrMLcXDR4IuIeSSeMo+hLw/cBJI0N359QQHg+Ij43jsrYAFl65IOs3ndar6th1he0X/z6k1OLZJ0paUPu/UhEjIzjkL8t6SGywTcfS03s9Ybvn9lsJ0UCwtWSrgDuBF7qFYuIB9quspmZ5e2JiAUT3McDwKsi4meSziNr1ZnHOIbvFwkIvwVcCLydXzUZRXpvZmY9FBH7cq/XSvq8pJmMY+h/kYCwFHh1mkLChpibjarN/QfDKY3W3BURIekMstsJngV+Qhq+DzxNNnz/j5rtq0hAeAg4Ct8AZjZuXhjHxkvSLWT3gM2UNApcQVqTJo3a/APgQ5IOAL8ALoiIAOoO/W92rCIBYRbwqKT1HNqH0HLYqZmZTUxELGuR/rdkw1LrpdUdvt9IkYBwRdGd2eAru9nI9yKMj5uLrBOK3Kn8T5JeBcyLiP8n6XCyyw8zMxsgRdZD+M/A/wW+kDbNIRvWZGZmA6TI5HYfBs4C9gFExPdpsAaBDQdPdlctbi6yTikSEF7IDzmVNJk21yYwGwT9slqa2XgVCQj/JOnjwCskvQP4O7J5jczMbIAUCQiXk62BsBn4INkQpv9eZqXMBkmZ9yC4ucg6qWVAiIiDEXFdRLwnIv4gvfa185Arsx/h3GO3lrZvM2us4bBTSbdFxB9K2kydPoOIOLnUmpm1yYvkmE1Ms/sQPpqetwJ/mtsu4NOl1cjMzHqiYUCIiJ3p5Wsj4ol8mqQ3lForMzPrumZNRh8C/hh4taRNuaTpwL+UXTEzM+uuZk1GNwNfA/6CbKTRmOciYm+ptTKzljzCyDqt4SijiPhpRDweEcsi4onco3AwkLRI0mOStkm6vE66JH0upW+SdHra/npJG3OPfZIuS2mflPR0Lu28cXxu6wDfsdyap722flJkttNxkTQJuAZ4B9nKPeslrYmI/ALPi8mWeptHttbntcCZEfEYcGpuP08Dq3PlPhsRnymr7mZmw6jIjWnjdQawLSK2p6kvbgWW1ORZAtwYmfuAoyTNrslzDvCD2o5tG2y+F8Gs+8oMCHOAp3LvR9O2dvNcANxSs21FamJaJenoegeXtFzSBkkb9u49WC+LmZnllBkQVGdb7Q1uTfNImgq8m2z+pDHXAq8ha1LaCVxV7+ARMRIRCyJiwYwZZX7M4eZ+hN5wh7KVocxvylFgbu79ccCONvMsBh6IiF1jGyJiV0S8GBEHgevImqbMzAZSagnZLenhBulLUovJxtQq8pZc2uOSNo+ltTpWmQFhPTBP0onpl/4FwJqaPGuA96fRRm8Gfpq7IQ5gGTXNRTV9DEuBun8kM7MBcT2wqEn63cApEXEq8J+AL9aknx0Rp0bEglYHKm2UUUQckLQCWEe25OaqiNgi6dKUvpJs5tTzgG3Az4GLx8qnpTrfQTbDat6nJZ1K1rT0eJ10s0rwkFPrhIi4R9IJTdJ/lns7jQmsV1NaQACIiLVkX/r5bStzr4NsRbZ6ZX8OvLLO9gs7XE2rqHOP3cq6Z07qdTWAai2O4/6D6pv0QuFzZmZNU85IRIy0ezxJS8luIj4G+L1cUgB3SgrgC632XWpAsOGw9MgHWb3vtF5Xw6wf7SnSlNNKRKwGVkt6K/Ap4HdT0lkRsUPSMcBdkh6NiHsa7cfDb8zMBkT6sn+NpJnp/Y70vJvs5t6mg3AcEMzM+pik10pSen06MBV4VtI0SdPT9mnAO2kxCMdNRjZQvEiODRpJtwALyfobRoErgCnwUp/sfyQbrflL4BfAeyMiJM0ia0aC7Lv+5oi4o9mxHBCsI9yPcKiyRhi5Q3n4RMSyFul/Cfxlne3bgVPaOZabjKzSPKeRWfc4IJiZGeCAYNY33FxkZXNAsI7xRHdm/c0BwczMAAcE6wP91rHsOYysXzkgWEe52agc7j+wbnBAMDMzwAHBrKUqzXRqViYHBLOKc3ORdYsDgnXcMPcjuEPZ+pkDgpmZAQ4I1ifaGXp6/Nw9JdbEbHA5IFgphrnZqJPcf2Dd5IBgZmZAyQFB0iJJj0naJunyOumS9LmUvimt9jOW9rikzZI25hehljRD0l2Svp+ejy7zM9j4+SphYnx1YN1WWkCQNAm4BlgMzAeWSZpfk20xMC89lgPX1qSfHRGn1ixCfTlwd0TMA+5O783MbILKvEI4A9gWEdsjYj9wK7CkJs8S4MbI3AccJWl2i/0uAW5Ir28Azu9gna3DhukqoZNDTn11YL1QZkCYAzyVez+athXNE8Cdku6XtDyXZ1ZE7ARIz8fUO7ik5ZI2SNqwd+/BCXwMM7PekbRK0m5JDzdIf19qct8k6V5Jp+TSmjbb1yozIKjOtto5AJrlOSsiTidrVvqwpLe2c/CIGImIBRGxYMYM9533UqeuEvpt1tPx8tWB1bgeWNQk/YfA2yLiZOBTwAgUbrY/RJnflKPA3Nz744AdRfNExNjzbmA1WRMUwK6xZqX0vLvjNTczq4iIuAfY2yT93oj4cXp7H9n3KBRrtj/E5A7Ut5H1wDxJJwJPAxcAf1STZw2wQtKtwJnATyNip6RpwK9FxHPp9TuBP8+VuQi4Mj1/pcTPYEPOE9vZeEx6IYr2Kc3Mj6IERiJiZAKHvgT4Wnpdr0n+zGaFSwsIEXFA0gpgHTAJWBURWyRdmtJXAmuB84BtwM+Bi1PxWcBqSWN1vDki7khpVwK3SboEeBJ4T1mfwTpn6ZEPsnrfab2uRuW5uWjo7KkZRTluks4mCwhvGdtUJ1vTXzhlXiEQEWvJvvTz21bmXgfw4TrltgOn1G5Pac8C53S2pmZm/UvSycAXgcXpOxKKNdsfwr2t1jXdHILaj/MZ+erAxkPS8cDtwIUR8b1c0kvN9pKmkjXbr2m2r1KvEMyGhae9trJIugVYSNbfMApcAUyBl1pcPgG8Evh8amY/kEZY1m22b3YsBwTrK+ceu5V1z5zU62p0nK8OrJGIWNYi/QPABxqkvazZvhk3GVlXDdOdy2b9xgHBzMwABwQzM0scEMx6zP0HVhUOCNZ1g9aP4BFGNigcEKzvDMskd2bd5oBg1oDnMbJh44Bg1kPuP7AqcUCwnuhGP0I/Tl9h1ksOCGZmBjggmE2IRxjZIHFAsJ6ZSLPRIIw0cv+BVY0DgpmZAQ4IZmaWOCBYT1X1ruWy70Fwc5FVkQOCmZkBDghWAWVeJZR5L8J4Rxj56sCqqtSAIGmRpMckbZN0eZ10SfpcSt8k6fS0fa6kb0jaKmmLpI/mynxS0tOSNqbHeWV+BquuQRhpZFYlpQUESZOAa4DFwHxgmaT5NdkWA/PSYzlwbdp+APiTiDgJeDPw4Zqyn42IU9Oj8PJwZmb9RtIqSbslPdwgve4P65T2uKTN6cfzhlbHKvMK4QxgW0Rsj4j9wK3Akpo8S4AbI3MfcJSk2RGxMyIeAIiI54CtwJwS62o9VtXO5U5zc5GNw/XAoibpjX5Yjzk7/Xhe0OpAZQaEOcBTufejvPxLvWUeSScApwHfyW1ekSLhKklH1zu4pOWSNkjasHfvwXF+BDOz3oqIe4C9TbLU/WE9nmOVGRBUZ1vtWL6meSQdAXwZuCwi9qXN1wKvAU4FdgJX1Tt4RIxExIKIWDBjhvvOrbM8ZYVVSLMf1gHcKel+Sctb7WhyCZXLV2pu7v1xwI6ieSRNIQsGN0XE7WMZImLX2GtJ1wH/2NlqW68sPfJBVu87ra0y5x67lXXPnFRSjTrLzUXDQ8/vL/rvPbOmbX8kIkbaPVydbWM/rM+KiB2SjgHukvRouuKoq8yfzuuBeZJOlDQVuABYU5NnDfD+1CnyZuCnEbFTkoAvAVsj4q/zBWouhZYCdTtarD8NS1+CWbJnrCUjPdoNBtDkh3VEjD3vBlaT9e02VFpAiIgDwApgHVmn8G0RsUXSpZIuTdnWAtuBbcB1wB+n7WcBFwJvrzO89NOp13wTcDbwX8v6DNYbgxgUfHVgJWr0w3qapOkAkqYB76TFD+gym4xIQ0LX1mxbmXsdwIfrlPsW9S+DiIgLO1xNq6DxNB91SqenrXAwsImQdAuwkKx5aRS4ApgCL32frgXOI/th/XPg4lR0FrA6a3BhMnBzRNzR7FilBgSziSgaFLrdj9BOh7KDgU1URCxrkd7oh/V24JR2juXhN1Zp/dx85GBg/cYBwSpvokGh22srT3101MHA+pIDgvWFfrlScCCwfuaAYNYG35Bmg8wBwQaCZz41mzgHBOsbVW82cnOR9TsHBDMzAxwQzMwscUCwvtKs2ajsfoRmHcpuLrJB4IBgfaeXQcFskDkgWF9qNygUvTmt0TxGvjqwYeCAYH2rW1cKDgY2LBwQrK+VHRQcDGyYOCBY3ysrKDgY2LBxQLCBUOSmtVb9CPn+AwcDG0YOCDYwyr6T2bOY2qBzQLCBUi8otNtsVO/qwIHAhoEDgg2cTl8pOBjYsHBAsIHUKCg06kdodP+Bg4ENk1IDgqRFkh6TtE3S5XXSJelzKX2TpNNblZU0Q9Jdkr6fno8u8zNY/xrPlUK+ucjBwKqgwPfo0ZJWp+/Q70p6U9GytUoLCJImAdcAi4H5wDJJ82uyLQbmpcdy4NoCZS8H7o6IecDd6b1ZXWNBod1+BAcDq4KC36MfBzZGxMnA+4Gr2yh7iDKvEM4AtkXE9ojYD9wKLKnJswS4MTL3AUdJmt2i7BLghvT6BuD8Ej+DmVkvFfkenU/245iIeBQ4QdKsgmUPMbnTtc+ZAzyVez8KnFkgz5wWZWdFxE6AiNgp6Zh6B5e0nOyqA+CFNx6/4+HxfIgOmgl0d7X3+qpQjy7XYUd6/mqP61FXFeoA1ahHFeoA8PqJ7mDfgR+tu+OZz88skPUwSRty70ciYiT3vsj36EPA7wPfknQG8CrguIJlD1FmQFCdbbU9d43yFCnbVPqjjgBI2hARC9op32lVqENV6lGFOlSlHlWoQ1XqUYU6jNVjovuIiEWdqAvFvguvBK6WtBHYDDwIHChY9hBlBoRRYG7u/XH86qdaqzxTm5TdJWl2ujqYDezuaK3NzKqj5fdoROwDLoZsoA7ww/Q4vFXZWmX2IawH5kk6UdJU4AJgTU2eNcD702ijNwM/Tc1BzcquAS5Kry8CvlLiZzAz66WW36OSjkppAB8A7klBosh38CFKu0KIiAOSVgDrgEnAqojYIunSlL4SWAucB2wDfk6Kco3Kpl1fCdwm6RLgSeA9Baoz0jpL6apQB6hGPapQB6hGPapQB6hGPapQB6hOPYp+j54E3CjpReAR4JJmZZsdTxFtNc2bmdmA8p3KZmYGOCCYmVnS1wGhjKkxSqrH+9LxN0m6V9IpubTHJW2WtHEiw90K1GGhpJ+m42yU9ImiZTtcjz/N1eFhSS9KmpHSOvW3WCVpt6S6955047woUIfSz4mC9Sj9vChQh26cE3MlfUPSVklbJH20Tp6ufF9UWkT05YOsk+QHwKvJhqk+BMyvyXMe8DWy8bhvBr5TtGyH6/E7wNHp9eKxeqT3jwMzu/C3WAj843jKdrIeNfnfBXy9k3+LtJ+3AqcDDzdI78Z50aoOpZ4TbdSjG+dF0zp06ZyYDZyeXk8HvteL74uqP/r5CqGsqTE6Xo+IuDcifpze3kc2HriTJvJ5uvq3qLEMuGWcx2ooIu4B9jbJUvp50aoOXTgnCtWjia79LWqUdU7sjIgH0uvngK1kd/LmdeP7otL6OSA0mvaiSJ4iZTtZj7xLyH6FjAngTkn3K5tuo8w6/LakhyR9TdIb2yzbyXog6XBgEfDl3OZO/C2K6MZ50Y4yzol2lH1eFNKtc0LSCcBpwHdqkqp2XnRdmXcql62nU2O0WY8so3Q22X/+t+Q2nxURO5TNyXSXpEfTL6pO1+EB4FUR8TNJ5wF/TzbLbE/+FmRNA/8SEflfjp34WxTRjfOiWEXKOyeK6sZ5UVTp54SkI8gCzmWR3bx1SHKdIj05L3qln68QJjI1RpGynawHkk4GvggsiYhnx7ZHxI70vBtYTXZ52vE6RMS+iPhZer0WmCJpZtH6d6oeORdQ0zTQob9FEd04L1oq+ZwopEvnRVGlnhOSppAFg5si4vY6WSpxXvRUrzsxxvsgu7rZDpzIrzp63liT5/c4tJPou0XLdrgex5Pdjf07NdunAdNzr+8FFpVUh2P51Y2IZ5Dd5a1u/y1Svt8ga1Oe1um/RW5/J9C4I7X086JAHUo9J9qoR+nnRas6dOOcSJ/pRuBvmuTpynlR5UffNhlFeVNjlFGPTwCvBD4vCeBAZLM6zgJWp22TgZsj4o6S6vAHwIckHQB+AVwQ2dne7b8FwFLgzoj411zxjvwtACTdQjZ6ZqakUeAKYEquDqWfFwXqUOo50UY9Sj8vCtQBSj4ngLOAC4HNymYFhWxhmeNz9Sj9vKg6T11hZmZAf/chmJlZBzkgmJkZ4IBgZmaJA4KZmQEOCGZmljggmJkZ4IBgZmaJA4INFUn/Ls11f5ikaWlu/Df1ul5mVeAb02zoSPqfwGHAK4DRiPiLHlfJrBIcEGzoSJoKrAeeJ5tL6MUeV8msEtxkZMNoBnAE2cpZh/W4LmaV4SsEGzqS1pCtenUiMDsiVvS4SmaV0LeznZqNh6T3k80serOkScC9kt4eEV/vdd3Mes1XCGZmBrgPwczMEgcEMzMDHBDMzCxxQDAzM8ABwczMEgcEMzMDHBDMzCz5/2tCe7OPX2O/AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "XX, TT = np.meshgrid(x, t)\n", "\n", "fig, ax = plt.subplots()\n", "CS = ax.contourf(XX, TT, T)\n", "ax.set_title('T(x,t)')\n", "ax.set_xlabel(\"x\") \n", "ax.set_ylabel(\"tiempo\")\n", "plt.colorbar(CS)" ] }, { "cell_type": "markdown", "metadata": { "id": "_0qpOOkTHUzk" }, "source": [ "Con esto hacemos una animación de la solución:\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "65gUIxPAK0sv", "outputId": "9cbd4a89-d8c1-4b11-b437-0fd1868ad36b" }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAApJElEQVR4nO3deXhU5d3/8fc3yUDCGiGAEpaIYFyQRdK6YbW1bXBpQWtxaW3tU2vV2u1po8XHaqut2tK6VGsRrbbWulUopXXB9qcVVFDCLgLKTggIBBIgmZDt/v1xBhzCTBbIzJnl87quuZKcZeabycl85pz7nvs25xwiIpK+MvwuQERE/KUgEBFJcwoCEZE0pyAQEUlzCgIRkTSnIBARSXMKAhGJCTO7xcweO8L7eNnMvt5RNUlkCoIkYWZ7w25NZhYM+/krftd3OMxsvZl91u869jMzZ2ZDW1h/tZm9GaPH7mVmfzezajPbYGZXtrDtcDObZWY7zOyQDwKZWYGZvWRmu8xsq5k9ZGZZYesnmtkKM9tjZu+b2YRY/E7Oubucc9eE1eTC64hQ98/M7Klm93G+c+7PsahPPqYgSBLOuW77b8BG4Athy/7qd33NtfQPn0yPEUe/B+qAfsBXgD+Y2clRtq0Hnge+GWX9w8A24BhgFHAOcAOAmeUDTwH/C/QASoCnzaxvh/wWkpycc7ol2Q1YD3w29H0G8BNgDVCB9wLRK7SuAHDAN4BNwC7gOuATwFKgEngo7H6vBt4CHgSqgJXAeWHrewJ/BLYAm4FfAJnN9r0P2BladxzwWqiuHcBfgdzQ9n8BmoAgsBe4CTgXKGvhd/0Z8ALeC9lu4JqWaorwvH0SmBv6vbcADwGdQutmh56r6lA9lzXb90SgFmgMra/swL9nV7wQOD5s2V+Ae1rZb6j3L3zI8hXABWE/TwYeCX1/GrCt2fbbgTNaO9bC/gZPNTu+vo735mQH8H9Rtt0Y2nZv6HZGs8cZF3oO6kPrl4SW/xe4JsIxVgmsBc4MLd+EF35fD7vPzsBvQo/9ETAFyPH7/zcRbzojSH7fAybgvevrj/di//tm25wGDAMuA+4H/g/4LHAyMNHMzmm27VogD7gdmG5mvULr/gw04L0AjQY+j/di3HzfvsAvAQPuDtV1IjAQ78UB59xVHHxm8+s2/r7j8cIgFy9YWqspXCPww9DvdgZwHqF3ys65T4W2GRmq57nwHZ1zK/BCdG5ofW6kBzCzh82sMsptaZS6jgcanXMfhC1bgvf3ORwPAJebWZfQGcD5wCuhdaXACjP7opllhi4L7cN7Y3C4xgKFeM/nbWZ2YoRt9j+/uaHnb274SufcK8BdwHOh9SOjPNZpoVp7A08Dz+K9sRkKfBV4yMy6hbb9Fd5zOyq0Ph+47bB+wxSnIEh+38Z7F1bmnNuH90J7abPLJnc652qdc6/iveN9xjm3zTm3GZiD9wK63zbgfudcfejFcBVwoZn1w3tB+YFzrto5tw3vndnlYfuWO+cedM41OOeCzrnVzrl/O+f2Oee2A/fiBdaRmOucm+Gca8K7tNFaTQc45xY45+aF6lsPPNIB9TR/jBucc7lRbiOi7NYN7wwsXBXQ/TDLeAMvRHYDZXgv/jNC9TUCT+K9iO4Lff22c676MB8L4Oehv/cSvACL9iLeEdY5554I/R7P4b25uCN0jL2Kd1Yx1MwM+BbwQ+fcTufcHrygiXhspLtUusaargYDfzezprBljXjXmvf7KOz7YISfu4X9vNk5F94AuQHvHf1gIABs8f7HAO+NxKawbcO/J3Td+XfA2Xgvahl4ZyxHIvwx2lJTeD3H44VREdAF7/hfcIT1dIS9eKEWrgewp713ZGYZwCy8kDsT72/7ON6745tCjfO/xrsMtxAYA8w0s/Odc4sPs/6tYd/XcPDx1NGaH7s45yIdz33w/sYLwo4NAzJjWFvS0hlB8tsEnN/snWd26N3+4ci3sP8cYBBQHnqcfUBe2OP0cM6FX75o3oPl7tCyEc65Hnin7tbC9tV4/7wAmFkm3j90uPB92lJTuD/gtXsMC9VzS7N6WtPqUL1mNqVZD6/w2/Iou30AZJnZsLBlI4Fo27ekF9675IdC75IrgCeAC0LrRwGznXOlzrkm59x84B28S4WRHPQ3AY4+jJqgDc9dG7dpqx14oXBy2LHR03mdLaQZBUHymwL80swGA5hZHzMbfwT31xf4npkFzOzLeNf2X3LObQFeBX5rZj3MLMPMjmvWvtBcd0INq6Fr1SXN1n8EDAn7+QMg28wuNLMAcCteg19Eh1FTd7zLJXvN7ATg+lbqae4jYICZdWqhputcWA+vZreIARW6LDMduMPMuprZWXhtIX+JtL15soFOoZ+zzaxz6L52AOuA680sy8xy8Rpzl4R2nw+cbWajQvuOxjtji9ZGsBivvSFgZkXApdF+91Zsx+sc0NrzWxA6qzkioUuHjwL37e8RZWb5ZlZ8pPedihQEye8BYCbwqpntAebhNagdrnfwGpZ34DX4Xhp6VwnwNbwXn/fxLvG8gNdFMZqfA6fiXe9+Ee/FLtzdwK2hhtQfO+eq8BpvH8PrAVSNd427Je2p6cfAlXiXXB7Fu8Yc7mfAn0P1TIyw/2t479K3mtmOVupqrxuAHLw2mmeA651zywHMbFDojGJQaNvBeO92958xBPHacva7BK8XznZgNV5j+g8BnHNvEOp9FTpepgF3ha6vR/JTvN5fu/D+nk8fzi/nnKvBO57eCj2/p0fY7G+hrxVmtvBwHqeZm/F+/3lmthv4D16jtjRjB18OlnRmZlfjddUb63ctIhI/OiMQEUlzCgIRkTSnS0MiImlOZwQiImku6T5QlpeX5woKCvwuQ0QkqSxYsGCHc67553KAJAyCgoICSktL/S5DRCSpmNmGaOt0aUhEJM0pCERE0pyCQEQkzSkIRETSnIJARCTNxazXkJkNxJsA42i8UQenOuceaLaN4Q2adgHeOOZXO+c6YrCpg8xYtJnJs1ZRXhmkf24OJcWFTBid39EPIyKSlGLZfbQB+JFzbqGZdcebIOLfzrn3w7Y5H2+ky2F4I2b+gSMbOfMQMxZtZtL0ZQTrGwHYXBlk0vRlAG0OAwWJiKSymAVBaKz4LaHv95jZCrw5Q8ODYDzwZGhGrHlmlmtmx4T27RCTZ606EAL7Besb+dnM5TgcOYFMOgcyyQlkkn3ga8aB5f9evpWf/uM9gvXeBGAKEhFJNXH5QJmZFeDNi/tOs1X5HDytYFlo2UFBYGbXAtcCDBo0iPYorwxGXF4ZrOeHzy2JuK41wfpGbp62lDc+2E7PnAA9cwLkdvFu3s+dDnz/xqpt3Dpj+RGdkYiIxFLMg8DMuuFNfvED59zu5qsj7HLIKHjOuanAVICioqJ2jZLXPzeHzRHC4Oge2Txz7enU1jcSrG+k9sCtiWBdI7UNjQTrGvnFiysi3u++hiZKN+yksqaePbUN7SnpwBnJ0T2zKejdlX49OnPw7JAH0xmFiMRSTIMgNN3gNOCvzrnms1OBdwYwMOznAXjz43aYkuLCg9oIAHICmfzk/BM4Nq9rq/s/8db6iEGSn5vDnJs+A0Bjk2N3sJ7KYD1VwXoqa+qoCn1/2z8iTztbGazn8qnzAMgOZFDQuyuDe3ehIK/rx9/37so7ayq4ZcZ7OqMQkZiJZa8hA/4IrHDO3Rtls5nAjWb2LF4jcVVHtg/Axy+Wh/uOOlqQlBR/PONdZoZxVNdOHNX10KlsH3ljbcQg6dejM7/58kjW76hmfUUNGyqqWb1tL6+v3E5dY1OLNQXrG5k8a5WCQEQ6RMzmIzCzscAcYBle91GAW4BBAM65KaGweAhvftUa4BvOuRZHlCsqKnLxHnTuSC7NNO+1BF6Q3H3JKRHvo7HJsaUqyPodNayvqObWGe9Fve8ffe54xg7LY8SAXDIzol9aEhExswXOuaKI65JtYho/guBIHUmQnHXPaxHPKAKZRn2j97frkZ3FmcflMXZYHmcPy2Nw79YveYlIelEQJLGWzijOHpbHW2sqePPD7bz54Q7Kq2oBGNgrh7FD8xg7tA9nDe3Nf1dtV2OzSJpTECS5tpxROOdYu6OaNz/cwZwPdzBvbQV793m9mcwg/M/c0qUpEUlNCoI01NDYxJKySq5+Yn7E7q39e2bz9qTzfKhMRPzQUhBo0LkUlZWZwZjBvdgb5TMO5VW1TJ29hj219XGuTEQSjYIgxfXPzYm4vHNWBne9tJIz736Ne15eybbdtXGuTEQShYIgxZUUF5ITyDxoWU4gk199aQQzbzyLTxX2YersNYz91evc/MJSVm/b61OlIuKXpJu8XtqntQ/U/f7KU9lQUc1jc9bxfOkmnivdxOdO6sd15wxhzOBefpYuInGixmI5oGLvPv48dwNPzl1PZU09RYOP4tvnHMd5J/Rl5pJydUEVSWLqNSTtUlPXwPPzN/HonHVsrgzSt3tndtXUHfgAG6gLqkiyUa8haZcunbK4+qxjeaPkXB64fBQ7qw8OAfh4vCMRSX4KAokqKzOD8aPyaWyKfNYYba4HEUkuCgJpVbQuqHndO8e5EhGJBQWBtCpSF1TwGpenzl5DU5QzBhFJDgoCadWE0fncfckp5OfmYHiT8vxywnA+d1I/7nppJVc9/g5bq/SBNJFkpV5Dcticczw3fxM//+f7dA5kcM8lIxg3/Gi/yxKRCNRrSGLCzLj8k4N48XtjGXhUF657agGTpi+lpq59cziLiL8UBHLEhvTpxrTrz+T6c4/j2fmbuOh3b7KsrMrvskSkjRQE0iE6ZWVw87gT+Os1p1FT18jFD7/FH/67JmrXUxFJHGojkA5XWVPHLX9fxkvLtnL6kF7cd9ko3lm7U0NUiPhIQ0xI3Dnn+NuCMn42cznOORqbHHUaokLEN2oslrgzMyYWDeSl751NQ7MQAA1RIZJIFAQSUwV5XWlo1BAVIolMQSAxF22IimjLRSS+FAQSc5GGqMgw+P55w3yqSETCKQgk5poPUXFUlwBNDmYs3kywrtHv8kTSnnoNiS+mLyzjR39bwhlDevPHr3+CnE6HDmonIh1HvYYk4Vxy6gDunTiSeWsr+J8/zdewFCI+UhCIby4ePYB7J47inXUKAxE/KQjEVxNG53PfZaN4d91OhYGITxQE4rvxoz4Og288oTAQiTcFgSSE/WEwf/1Orn5iPtX7FAYi8RKzIDCzx81sm5m9F2V9TzP7p5ktMbPlZvaNWNUiyWH8qHzuv3w0peu9MwOFgUh8xPKM4E/AuBbWfwd43zk3EjgX+K2ZdYphPZIEvjiyPw9cPpoFG3cpDETiJGZB4JybDexsaROgu5kZ0C20rf7rhS+M7M8Dl49iwcZdXP3Eu+xVGIjEVJaPj/0QMBMoB7oDlznnmiJtaGbXAtcCDBo0KG4Fin8uGtEfw/jes4u4+vF3uXTMAB58bbXmMxCJAT+DoBhYDHwGOA74t5nNcc7tbr6hc24qMBW8TxbHs0jxz4UjjgHgxqcXsnDjLvZPdra5Msik6csAFAYiHcDPXkPfAKY7z2pgHXCCj/VIArpwxDEHxiYKp/kMRDqOn0GwETgPwMz6AYXAWh/rkQS1q6Y+4nLNZyDSMWJ2acjMnsHrDZRnZmXA7UAAwDk3BbgT+JOZLQMMuNk5tyNW9Ujy6p+bw+YIL/qaz0CkY8QsCJxzV7Syvhz4fKweX1JHSXEhk6YvI1j/8ZDVnTIzKCku9LEqkdThZ2OxSJvsbxCePGsV5ZVBMjOMzAwYM/gonysTSQ2aj0CSzvod1XzhoTcZ1KsL064/k+yA5jIQaY3mI5CUUpDXlQcuH8X7W3Zzy/RlJNubGZFEoyCQpPSZE/rxg/OOZ/qizTw5d4Pf5YgkNQWBJK3vfmYonz2xL3f+633mr29pNBMRaYmCQJJWRoZx72WjGNirCzf8dSEf7a71uySRpKQgkKTWIzvAI1eNoXpfA9c/tYC6hojDVYlICxQEkvSO79edyZeOZOHGSu7413K/yxFJOgoCSQkXjjiGb58zhKfmbeT50k1+lyOSVBQEkjJKPl/IWUN7c+uM91haVul3OSJJQ0EgKSMrM4MHrziVPt06c91fFlCxd5/fJYkkBQWBpJReXTvxyFVjqKiu47vPLKKhUY3HIq1REEjKGZ7fk19efApvr6ngV6+s9LsckYSnQeckJV06ZgBLyyp5dM46RgzI5Qsj+/tdkkjC0hmBpKxbLzyJosFHcdMLS1m59ZAZUEUkRKOPSkrbtruWix58k8amJjplZbK1qpb+uTmUFBdqvmNJKxp9VNJW3x7ZXPHJgVRU17OlqhYHbK4MMmn6MmYs2ux3eSIJQUEgKe+FBYe+4AfrG5k8a5UP1YgkHgWBpLxok9xHWy6SbhQEkvKiTXIfbblIulEQSMorKS4kp9l0llkZRklxoU8ViSQWfY5AUt7+3kGTZ62ivDJITqdMauoaGdiri8+ViSQGdR+VtLOntp5x988hkGm89P2z6dJJ74ck9an7qEiY7tkBfjtxJBt21nDPyxqCQkRBIGnp9CG9+Z+zjuXJuRuY/cF2v8sR8ZWCQNJWSXEhQ/t246YXllJVU+93OSK+URBI2soOZHLfxFHs2LuP22e+53c5Ir5REEhaO2VAT777mWHMWFzOi0u3+F2OiC8UBJL2bvj0cYwc0JNbZyxj255av8sRiTsFgaS9QGYGv504ipq6RiZNW0aydakWOVIxCwIze9zMtplZ1IuvZnaumS02s+Vm9kasahFpzdC+3bh53An8v5XbeL50k9/liMRVLM8I/gSMi7bSzHKBh4EvOudOBr4cw1pEWnX1mQWcMaQ3d/zzfTbtrPG7HJG4iVkQOOdmAztb2ORKYLpzbmNo+22xqkWkLTIyjN9MHEmGGT96fgmNTbpEJOnBzzaC44GjzOy/ZrbAzL4WbUMzu9bMSs2sdPt2ffhHYic/N4fbv3gy767fyeNvrvO7HJG48DMIsoAxwIVAMfBTMzs+0obOuanOuSLnXFGfPn3iWaOkoS+dms/nT+rH5Fmr+OCjPX6XIxJzfgZBGfCKc67aObcDmA2M9LEeEQDMjLsuOYXu2Vn88LnF1DU0+V2SSEz5GQT/AM42sywz6wKcBqzwsR6RA/K6deaXF5/C8vLdPPTah36XIxJTMRt/18yeAc4F8sysDLgdCAA456Y451aY2SvAUqAJeMw5p8/5S8IYN/xovnTqAH7/3zV8+oS+jB50lN8licSE5iMQacHu2nrG3Teb+sYmApkZbKmqpX9uDiXFhQcmvBFJBpqPQOQw9cgOMGF0Ptv31lFeVYsDNlcGmTR9GTMWbfa7PJEOoSAQacU/FpcfsixY38jkWat8qEak4ykIRFpRXhls13KRZKMgEGlF/9ycdi0XSTYKApFWlBQXkhPIPGhZp6wMSooLfapIpGPFrPuoSKrY3zto8qxVlFcGyTCjZ3YW44Yf7XNlIh1DQSDSBhNG5x8IhDkfbueqP77L/f/5kJ+cf4LPlYkcOV0aEmmns4f14bKigUydvYalZZV+lyNyxBQEIofhlgtPpE/3ztz0wlKNRSRJT0Egchh65gS46+JTWLl1Dw//d7Xf5YgcEQWByGE678R+jB/Vn4deW83Krbv9LkfksEUNAjNTQ7JIK27/wsn0zAlw0wtLaWjUJSJJTi2dEbwbtypEklSvrp24Y/xwlpZV8ZhmNJMk1VIQWNyqEEliF5xyNMUn9+Pef3/Amu17/S5HpN1auvzTx8z+N9pK59y9MahHJOmYGXeOH87n7pvNzS8s5flvn0FGht5HSfJo6YwgE+gGdI9yE5GQvj2yue2ikyjdsIsn5673uxyRdmnpjGCLc+6OuFUikuQuOTWffy4t51evrOK8E/sxsFcXv0sSaRO1EYh0EDPjrotPITPDuHnaUpJt9j9JXy0FwXlxq0IkRfTPzWHSBSfw9poKnp2/ye9yRNokahA453bGsxCRVHHFJwZxxpDe3PXiCrZUafIaSXz6ZLFIB8vIMO750ik0NDlumb5Ml4gk4SkIRGJgcO+u/Li4kNdXbWfGYk1yL4lNQSASI1efWcCpg3L5+T/fZ/uefX6XIxKVxhMSiZHMDOPXl47kgt/N4Zo/v8uOvfWUVwbpn5tDSXHhgYluRPymMwKRGBratxufP6kvS8p2s7kyiAM2VwaZNH0ZMxbpkpEkBgWBSIwt3Fh5yLJgfSOTZ62KfzEiESgIRGJsS2VtxOXllepaKolBQSASY/1zc9q1XCTeFAQiMVZSXEhOIPOgZdlZGZQUF/pUkcjB1GtIJMb29w6aPGsV5aEG44K8rowf1d/fwkRCFAQicTBhdP6BQPjz2+u5feZynn53I185bbDPlYnE8NKQmT1uZtvM7L1WtvuEmTWa2aWxqkUkkVx1+mDGDs3jly+uYENFtd/liMS0jeBPwLiWNjCzTOBXwKwY1iGSUDIyjF9fOoLMDONHzy+hsUljEYm/YhYEzrnZQGsjmH4XmAZsi1UdIomof24Od4w/mdINu3h0zlq/y5E051uvITPLBy4GprRh22vNrNTMSrdv3x774kTiYMKofM4ffjT3vvoBK7bs9rscSWN+dh+9H7jZOdfY2obOuanOuSLnXFGfPn1iX5lIHJgZv5gwnB45WfzwucXsa2j1X0EkJvwMgiLgWTNbD1wKPGxmE3ysRyTuenfrzN2XjGDl1j088J8P/S5H0pRvQeCcO9Y5V+CcKwBeAG5wzs3wqx4Rv3zupH5MLBrAlDfWsGDDLr/LkTQUy+6jzwBzgUIzKzOzb5rZdWZ2XaweUyRZ/fSikzimZw4/en4xNXUNfpcjaSZmHyhzzl3Rjm2vjlUdIsmge3aA304cyRWPzuPul1Zy54ThfpckaURjDYkkiNOH9OZ/zjqWv8zbwBsfqHecxI+CQCSBlBQXMrRvN256YQlVNfV+lyNpQkEgkkCyA5ncN3EUFXvruH1mi6OziHQYBYFIgjllQE+++5lhzFhczotLt/hdjqQBBYFIArrh08cxYkBPbp2xjG27I89wJtJRFAQiCSiQmcG9E0dSU9fIT6YvwzkNTCexoyAQSVBD+3bn5nEn8NrKbTw3f5Pf5UgK08Q0Igns6jML+Pf7H3HbP97jvv98wLbd++ifm0NJceGBiW5EjpTOCEQSWEaG8fmT+1HX6Pho9z4csLkyyKTpy5ixaLPf5UmKUBCIJLjH5qw7ZFmwvpHJs1b5UI2kIgWBSIIrrwy2a7lIeykIRBJc/9ycdi0XaS8FgUiCKykuJCeQedAyA75/3lB/CpKUoyAQSXATRudz9yWnkJ+bgwG9u3bCAW+urtDnC6RDqPuoSBKYMDr/oO6iv399NZNnrWLEgJ5cc/YQHyuTVKAzApEkdMO5x1F8cj/ufnklc9dU+F2OJDkFgUgSMjN+8+WRFPTuwo1PL1QPIjkiCgKRJNU9O8AjVxWxr6GJ659aQG19o98lSZJSEIgksaF9u/GbL49kSVkVP5u53O9yJEkpCESS3LjhR/OdTx/Hs/M38fQ7G/0uR5KQgkAkBfzv5wr51PF9uH3meyzcuMvvciTJKAhEUkBmhvG7y0dxdM9sbnhqIdv37PO7JEkiCgKRFJHbpROPfLWIymAd3/nrQuobm/wuSZKEgkAkhZzUvwf3XDKCd9fv5K6XVvhdjiQJfbJYJMVMGJ3PkrJKnnhrPSMG9OTi0QP8LkkSnM4IRFLQLRecyGnH9mLS9GUsL6/yuxxJcAoCkRQUyMzgoStPJTenE9c9tYDKmjq/S5IEpktDIimqT/fO/OGrp3LZI/O47JG57NnXwJbKWs15LIfQGYFIChs96CgmjO7Pqo/2Ul5ZqzmPJSIFgUiKe2v1oaOTas5jCRezIDCzx81sm5m9F2X9V8xsaej2tpmNjFUtIulMcx5La2J5RvAnYFwL69cB5zjnRgB3AlNjWItI2tKcx9KamAWBc242sLOF9W875/YPijIPUGdnkRiINOcxwLmFeT5UI4koUdoIvgm8HG2lmV1rZqVmVrp9+/Y4liWS/JrPedy/ZzbD+nbl6Xc3MW1Bmd/lSQLwvfuomX0aLwjGRtvGOTeV0KWjoqIizdYt0k7N5zwO1jXyrSdL+fELS3DApWN0Qp7OfD0jMLMRwGPAeOecJl4ViZOcTpk89vUizjouj5IXlvC30k1+lyQ+8i0IzGwQMB24yjn3gV91iKSr7IAXBmOH5nHTtKU8rzBIWzG7NGRmzwDnAnlmVgbcDgQAnHNTgNuA3sDDZgbQ4JwrilU9InKo7EAmj36tiG89WcrN05aCg4mfGOh3WRJnMQsC59wVray/BrgmVo8vIm2zPwyu/csCbp6+FFAYpJtE6TUkIj7KDmQy9aoxnD2sDzdNW8pz8zX3cTpREIgI8HEYfOr4Ptw8bZnCII0oCETkgP1hcE4oDJ59V2GQDhQEInKQ7EAmj1w1hnML+/CT6ct4RmGQ8nz/QJmIJJ7sQCZTvjqG659awKTpy1i4cRdvr66gvDKo+QxSkM4IRCSi7EAmU64aw4lHd+dvpWVsrgxqPoMUpSAQkag6Z2VSFaw/ZLnmM0gtCgIRadGWqtqIyzWfQepQEIhIi6LNW5DbJYBzGgMyFSgIRKRFkeYzyDDYVVPP959dHPHSkSQX9RoSkRbt7x00edaqA72GfvS54ymvCnLffz5kwYZd3H/5KD5R0MvnSuVwWbKd2hUVFbnS0lK/yxARYNHGXfzgucVs2lnDdz49lO+dN4xApi40JCIzWxBtYE/9xUTksI0edBQvfu9sLjl1AA++tpovT5nLhopqv8uSdlIQiMgR6dY5i998eSQPXTmatdv3csEDc/hb6SY1JCcRBYGIdIiLRvTnlR98iuH5PSl5YSk3PrOIqho1JCcDNRaLSIfpn5vD0986nUdmr+HeVz9g0YZd3HvZKLZW1R7U2KwhKhKLGotFJCaWbKrkB88tZt2OarIyjIamj19rcgKZ3H3JKQqDOFJjsYjE3ciBufzru2Pp0inzoBAADVGRaBQEIhIzXTtnEaxrjLhOQ1QkDgWBiMRUtCEqMjOMvy8qo76xKc4VSXMKAhGJqUhDVAQyjbxunfjhc0s459ev88c311G9r8GnCkW9hkQkpiINUVFSXMgXR/bnvx9sY8oba7nzX+/zu//3IVedPpivn1lAn+6dfa46vajXkIj4buHGXUx9Yy2z3t9KIDODS8cM4FtnD+HYvK5+l5YyWuo1pCAQkYSxdvteHp2zjmkLvbaDcScfzbfPOY71O6r1OYQjpCAQkaSybU8tf357PX+Zu4HdtQ1kGIT3QNXnENpPnyMQkaTSt3s2JcUn8Pak8+iZk0WzjyEQrG/k7pdX+FNcClJjsYgkrG6ds9gdjNyb6KPd+zh38uuMHZbH2KF9OOO43vTMCcS5wtSgIBCRhNY/N4fNET581jMniyF9ujF94WaemreRDPM+zXz20DzGDuvD6EG5B+ZGmLFos9oYWqA2AhFJaDMWbWbS9GUE6z/+hHJ4G0FdQxOLNu7izdU7eHP1DpZsqqTJQddOmZw+pDc9srN46b2t7Gtoirh/ulBjsYgktfa8o68K1jN3TQVvrt7Omx/uYH1FTcTtenftxPQbziQ/N4esVmZVS4UzCgWBiKStY3/yIi29ymVlGAN7dWFw7y4U9O5KQe8uDM7ryrG9u5J/VA4vLt3S4hlJWyRCkLQUBDFrIzCzx4GLgG3OueER1hvwAHABUANc7ZxbGKt6RCQ9RWtjyOvWiZuKT2B9RTUbKmpYt6Oa+et2Uh02SF5mhgHQGGH01Dv/9T7H5nUlt0uAnjkBumcHDmwfrvmlrc2VQSZNXwaQMEESy8biPwEPAU9GWX8+MCx0Ow34Q+iriEiHKSkujPiO/tYLTzrkxdQ5x469dWyoqGbdDi8gHnp9dcT7raiuY/zv3zrwsxn0yPZCYX849MwJ8PrKbQc9NnwcJH17dCY7kElOIDPsawbZgUw6Z2VgZh0SJK2JWRA452abWUELm4wHnnTetal5ZpZrZsc457bEqiYRST/RxjqK9CJqZvTp3pk+3TtTVNALgL8v2hz1jOJXXxpBZU09VcF6KoP17A7WU1lTR2XQW7Z5V/CgM4xwFdV1XPnoO1HrNoPsrEz2NTRG/BzF5FmrEj8I2iAf2BT2c1lo2SFBYGbXAtcCDBo0KC7FiUjqmDA6/7BfNFs6ozjvxH6t7n/WPa9FDJI+3Trz4JWjCdY3sq++kWB9I7X1TQTrGqltaKS2zlv26Jx1Ee+3I+dz8DMIDr2YRuQ2HefcVGAqeI3FsSxKRCRce84oIokWJP934YmcPqR3q/u/tGxrxCCJNs/D4fAzCMqAgWE/DwDKfapFRCSqIzmjiFWQlBQXHlY9kfgZBDOBG83sWbxG4iq1D4hIKvIzSNoilt1HnwHOBfLMrAy4HQgAOOemAC/hdR1djdd99BuxqkVEJJkdSZC0RSx7DV3RynoHfCdWjy8iIm2jYahFRNKcgkBEJM0pCERE0pyCQEQkzSXd6KNmth3YcJi75wE7OrCcjpKodUHi1qa62kd1tU8q1jXYOdcn0oqkC4IjYWal0YZh9VOi1gWJW5vqah/V1T7pVpcuDYmIpDkFgYhImku3IJjqdwFRJGpdkLi1qa72UV3tk1Z1pVUbgYiIHCrdzghERKQZBYGISJpLmSAws3FmtsrMVpvZTyKsNzP7XWj9UjM7ta37xriur4TqWWpmb5vZyLB1681smZktNrPSONd1rplVhR57sZnd1tZ9Y1xXSVhN75lZo5n1Cq2L5fP1uJltM7P3oqz36/hqrS6/jq/W6vLr+GqtrrgfX2Y20MxeN7MVZrbczL4fYZvYHl/OuaS/AZnAGmAI0AlYApzUbJsLgJfxZkY7HXinrfvGuK4zgaNC35+/v67Qz+uBPJ+er3OBfx3OvrGsq9n2XwBei/XzFbrvTwGnAu9FWR/346uNdcX9+GpjXXE/vtpSlx/HF3AMcGro++7AB/F+/UqVM4JPAqudc2udc3XAs8D4ZtuMB550nnlArpkd08Z9Y1aXc+5t59yu0I/z8GZqi7Uj+Z19fb6auQJ4poMeu0XOudnAzhY28eP4arUun46vtjxf0fj6fDUTl+PLObfFObcw9P0eYAXe/O3hYnp8pUoQ5AObwn4u49AnMto2bdk3lnWF+yZe6u/ngFfNbIGZXdtBNbWnrjPMbImZvWxmJ7dz31jWhZl1AcYB08IWx+r5ags/jq/2itfx1VbxPr7azK/jy8wKgNHAO81WxfT48nOqyo5kEZY17xcbbZu27Hu42nzfZvZpvH/UsWGLz3LOlZtZX+DfZrYy9I4mHnUtxBubZK+ZXQDMAIa1cd9Y1rXfF4C3nHPh7+5i9Xy1hR/HV5vF+fhqCz+Or/aI+/FlZt3wgucHzrndzVdH2KXDjq9UOSMoAwaG/TwAKG/jNm3ZN5Z1YWYjgMeA8c65iv3LnXPloa/bgL/jnQbGpS7n3G7n3N7Q9y8BATPLa8u+sawrzOU0O22P4fPVFn4cX23iw/HVKp+Or/aI6/FlZgG8EPirc256hE1ie3x1dMOHHze8M5u1wLF83GBycrNtLuTgxpZ327pvjOsahDdv85nNlncFuod9/zYwLo51Hc3HHzj8JLAx9Nz5+nyFtuuJd523azyer7DHKCB642fcj6821hX346uNdcX9+GpLXX4cX6Hf+0ng/ha2ienxlRKXhpxzDWZ2IzALrxX9cefccjO7LrR+CvASXsv7aqAG+EZL+8axrtuA3sDDZgbQ4LzRBfsBfw8tywKeds69Ese6LgWuN7MGIAhc7rwjz+/nC+Bi4FXnXHXY7jF7vgDM7Bm8ni55ZlYG3A4EwuqK+/HVxrrifny1sa64H19trAvif3ydBVwFLDOzxaFlt+CFeFyOLw0xISKS5lKljUBERA6TgkBEJM0pCERE0pyCQEQkzSkIRETSnIJARCTNKQhERNKcgkDkCJnZJ0JjxGebWdfQmPLD/a5LpK30gTKRDmBmvwCygRygzDl3t88libSZgkCkA5hZJ2A+UIs3rk+jzyWJtJkuDYl0jF5AN7wZprJ9rkWkXXRGINIBzGwm3uxQxwLHOOdu9LkkkTZLidFHRfxkZl/DG9XzaTPLBN42s884517zuzaRttAZgYhImlMbgYhImlMQiIikOQWBiEiaUxCIiKQ5BYGISJpTEIiIpDkFgYhImvv/dcurt9AWYlgAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def pintar_solucion(T_n, n):\n", " # Clear the current plot figure\n", " plt.clf()\n", " plt.title(f\"Temperature at t = {n*dt:.3f} unit time\")\n", " plt.xlabel(\"x\")\n", " plt.ylabel(\"T\")\n", " plt.plot(x,T_n,'-o')\n", " return plt\n", "\n", "def animate(n):\n", " pintar_solucion(T[n], n)\n", "\n", "anim = animation.FuncAnimation(plt.figure(), animate, interval=100, frames=nt, repeat=False)\n", "#anim.save(\"heat_equation_solution.gif\")\n", "\n", "HTML(anim.to_html5_video()) " ] }, { "cell_type": "markdown", "metadata": { "id": "-LX2QeQwJj-1" }, "source": [ "Ahora debemos comprobar el resultado que hemos obtenido. Para ello, podemos compararlo con la solución analítica para este mismo problema:\n", "\n", "$$ T(x,t)=1.5-0.5erf\\left(\\frac{x-1}{2\\sqrt{bt}}\\right)$$\n", "\n", "donde $erf()$ es la función de error que se calcula utilizando ```special.erf()```. Se pide:\n", "\n", "1.d) Calcula la solución exacta en el tiempo final, almacénala en un array ```T_exacta``` y represéntala gráficamente junto a la solución numérica.\n", "\n", "1.e) Calcula el error de la aproximación numérica en cada punto y almacénalo en un vector:\n", "\n", "$$ \\epsilon_i=\\left|T_i-T_{exacta}(x_i)\\right| $$\n", "\n", "y después calcula el máximo de ese vector\n", "\n", "$$L_{\\infty}= \\max(\\epsilon)$$\n", "\n", "A este valor lo llamamos \"norma infinito del error\" y nos da una medida del error máximo cometido por la aproximación numérica.\n", "\n", "1.f) Repite este cálculo para distintos valores de N y explica: ¿Cómo cambia el error al variar N? (no hay que escribir más codigo, simplemente prueba varios valores de N ejecutando lo que ya tienes y explica los resultados)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "f7Fkm69VciSH", "outputId": "5630d8d0-06d2-4e69-867a-1fb014b63a91" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.3\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAE9CAYAAAAGZmUpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABHS0lEQVR4nO3dd3hUVf7H8feZksykTSCNklClQ0JVRCmCKFJUdhV0FcGfZdFlAV1dbIuIay8o6ypiL9hQcIFFV6WIYqNI7yBCaOmN9Jnz+2NCIKSHmbmZ5Pt6njxJ7rnlm+Eyn7ntHKW1RgghRONlMroAIYQQxpIgEEKIRk6CQAghGjkJAiGEaOQkCIQQopGTIBBCiEbOYnQBtRUZGanbtGljdBlCCOFXNmzYkKK1jqqoze+CoE2bNqxfv97oMoQQwq8opX6vrE1ODQkhRCMnQSCEEI2cBIEQQjRyfneNQIjGqKioiMTERPLz840uRdRzNpuN2NhYrFZrjZeRIBDCDyQmJhIaGkqbNm1QShldjqintNakpqaSmJhI27Zta7ycnBoSwg/k5+cTEREhISCqpJQiIiKi1keOEgRC+AkJAVETddlPvBYESqk4pdQqpdROpdR2pdS0CuZRSqm5Sql9SqktSqneXilmyycwpzvMCnd/3/KJVzYjREN2/PhxrrvuOtq3b0/Xrl0ZOXIke/bsOef1rl69mtGjRwOwZMkSnnzyySrnHzBgwDlvE2DmzJl88803Vc7jqXpCQkJqVZuvefMaQTHwN631RqVUKLBBKfW11nrHGfNcAXQo+boAeKXku+ds+QSWToWiPPfvmYfdvwPEj6v5OlbMhsxEcMTCsJk1X1aIBkBrzdixY5k4cSIfffQRAJs2beLEiRN07NjRY9u58sorufLKK6uc54cffvDItmbPnl2v6jGS14JAa30MOFbyc7ZSaifQEjgzCK4C3tXuYdJ+UkqFK6WalyzrGStmnw6BU4ryyFs2gzXHg9GWYEb0bg8BwfyUWMBvGUVlZm1/fDnnb51VJkiKP/8rv+xP5feWowgOtHBlQgsAVu1K4nhW2XNz4XYrV+jvYMVsdGYiubZmbO40jd9bjgIgMiSQ4V1jAPjvlmNk5ZfdfjOHjUs6RQPwn01HyC10lmmPaxLExR0iAfh0QyJFTleZ9raRwfRvFwHAR78c4uzx6DrGhNCndVOKnS4Wbkgs9/J1bR5GQlw4+UVOFv96pFx7fKyDbi0cZOcXsWxL+X+2Pq2b0DEmlPSThXy5/Xi59gvaNqVdVAhJ2fms2JlUrv3i8yKJaxrE0Yw8vt2TXK59SKcomjvsHErNZe3+lHLtl3aJISo0kP3JOfzyW1q59pHdm+MIsrL7eDYbD6WXa78yoQXBgRa2Hclkz4lsRvZojs1qLjdfQ7dq1SqsViuTJ08undazZ0/AHRJ///vf+eKLL1BK8dBDDzF+/HhWr17NrFmziIyMZNu2bfTp04f3338fpRRffvkl06dPJzIykt69T58IePvtt1m/fj0vvfQSJ06cYPLkyRw4cACAV155hQEDBhASEkJOTk6dtnumSZMmMXr0aK655hratGnDxIkTWbp0KUVFRSxcuJDOnTvXqp6cnByuuuoq0tPTKSoq4p///CdXXXWVl/9lPMMndw0ppdoAvYCfz2pqCRw+4/fEkmll3lGUUrcDtwO0atWqdhvPdL+5rbbbeTiqKXaXxq5dBLk09gPTCXK5WL1TY3e5CNKaQBdYXGbM2ozJaSWVTNZaNEGmAIJcmkinkyaufLps+idvrUvCHBrFlXEjIDiKN747wPf7U8tsfnKTDVxR/DIU5aGA4Pxj9Nw0k4/WHWKJ62L6tG5SGgQvfLOHvUk5ZZYf1DGqNAie+mIXRzPLBs3IHs1Kg2D20u1k5ReXab+2T2xpEDz0+TaKXWWjYNKANu4gcGnuX7S13Mv3l0vakxAXzsmC4grb/z6iE91aOEg/WVRh++yrutExJpRjmfkVtj93bQLtokI4lJpbYfu8G/sQ1zSIPSeyK2x/75bzae6ws/VIZoXtne4MJSo0kA0H0yts79O6CY4gK2v3pTB72Y5y7YM6RhEcaGH17iSe/WoPARYTo+NblJuvoTv1hlqRRYsWsWnTJjZv3kxKSgr9+vVj0KBBAPz6669s376dFi1acNFFF7F27Vr69u3LbbfdxsqVKznvvPMYP358heudOnUqgwcPZvHixTidTnJyyv7fqO12L7744ir/xsjISDZu3MjLL7/Ms88+y+uvv16remw2G4sXLyYsLIyUlBT69+/PlVde6RfXdrweBEqpEOAzYLrWOuvs5goWKTeIstZ6PjAfoG/fvrUbZNkRC5mHiXIWc+nJXPKUiVyT4qQlkJymbUlx5lOgC8hzFpLrKiJPF5+1gqblVmnRmginkyjnAqKKnTy64GminE6u1iZuaRFCREA4Te2ROOwxBO9bXu6IJEgV8nzEEh649WGs5tMvwQe39cd51ht1gOX0ZZzPp1yEq+wHfgLPaP/67sGcPQS1/YxPr9/PGFrub7EHmEvX89P9w8q1BwW628ODAipsD7G5d6EW4bYK20NL2jvEhFTY7rC773XuEeuosD08yN3ev11Ele3DukRX2N40OACA0QnNGdSxfH9bESHu9vH94hjZo3m59siS9qt6tuTZr/aQdrKw3DxGGP/qj+WmjY5vzoQL25BX6GTSW7+Ua7+mTyzX9o0j7WQhd7y/oUzbx3++sM61fP/991x//fWYzWZiYmIYPHgw69atIywsjPPPP5/Y2FjAfQRx8OBBQkJCaNu2LR06dADgxhtvZP78+eXWu3LlSt59910AzGYzDofjnLZbXRD84Q9/AKBPnz4sWrSo1vVorXnggQdYs2YNJpOJI0eOcOLECZo1a1bta2g0rwaBUsqKOwQWaK3Lv7LuI4C4M36PBY56tIhhM2HpVLoV5tEtteTQ32qHMU9UeJ7fpV3kF+eTW5xLXnEeuW9dQd7JE+QqEzkmRYrZTLLFTLItlOSWvUjMS2FTQQbpzjPf7LOALMy5+4mICSXKGURUsZMYp5OOhUV0KSikY04izdbOhJiuEN0NorsQFVrJBaWSaxTR1VyjiAmzVflSNHNU3q6UqrLdbKq63WI2VdluraY90GKmmaPyUy4267m1BwVYCAqofHcPDrQQHFh5e3RYIABZeUWVztOQdevWjU8//bTCNn32p48zBAYGlv5sNpspLnZ/0PLEp+S6bLcqp5ap6fxnW7BgAcnJyWzYsAGr1UqbNm385gFArwWBcv9LvwHs1Fo/X8lsS4ApSqmPcF8kzvTo9QE4/YZZw4u9JmUiyBpEkDXIPWHIzLIXm8EdJJc8XGYdRc4iUvJSSMpLIiXX/T05N5mU9fNJchZxzGJhg83Gx2HuT/BmrWl/dCldDyykS0EhXQsL6RTUHHt0d4jpdjogjmyE/04/t4vd4pwFWszYrCYy60kQVPUJ3h5grrK9aXBArY8Ahg4dygMPPMBrr73GbbfdBsC6devIzc1l0KBBvPrqq0ycOJG0tDTWrFnDM888w65duypcV+fOnfntt9/Yv38/7du358MPP6xwvmHDhvHKK68wffp0nE4nJ0+eJCwsrLS9tts9V9XVk5mZSXR0NFarlVWrVvH775V29lnvePOI4CJgArBVKbWpZNoDQCsArfU8YDkwEtgH5AI3e6WS+HF1f9OsYZBYzVaahzSnechZpxcszUqDRANHLWZ22kPY0W0kO8ywJmUbnxdmAmBC0y5vG112/kTXTYV0KSikc2EhwWd/8inKc9cjQeBT8S3DS09lNTZKKRYvXsz06dN58sknsdlstGnThhdeeIFBgwbx448/kpCQgFKKp59+mmbNmlX6hmyz2Zg/fz6jRo0iMjKSiy++mG3btpWb78UXX+T222/njTfewGw288orr3DhhacDbOzYsbXa7rmqrp4bbriBMWPG0LdvX3r27Ennzp29Uoc3qKoOr+qjvn37ar8bj6CK20+11pzIPcGO1B3sTNvp/p66g+Q89x0wSmtaFxXTtbCQ7gWFDMvNpUVxyZ1D0zZDkzYG/VHCl3bu3EmXLl2MLkP4iYr2F6XUBq1134rml76GfKGKIxKlFM2Cm9EsuBlDW52+mJucm8zOtJ1sX3onO3URG2yBLA8J5umIJvTOz2d0zkku+1cvHM16Qrex0O1qCK/lHVVCCIEEQb0VFRRFVFAUgy46dY0ihUSLmS+Cg1kWGsLsyAgejzQxyJnJ6LWPM+ibfxDYsp87FLpe5T7yEB711Je7OJyWy0t/8s4D8EIYRYKgvjvjGkVsZiK34eDWfv9gV2w8yw4s44vfvmCluYhQUwDDC9MZ/e0j9PnfA5jiLoBuf3CHwsHv5MloDziemc+mwxlGlyGEx0kQ+IOzTi0poAvQJaILd/e5m1+O/8KyA8v48vdvWNQ8hhhzECPzUxm14iE6fTkDlAl0yQMIctdRnYXZLPXmriEhPEl6H/VzZpOZC1tcyGMXP8bq8at5ZtAzdG7ej/cCXFwT25w/xLbkzbBgjpvPuMf+1F1HolYcdis5BcW4XP51g4UQ1ZEjggbEbrEzou0IRrQdQXp+Ov87+D+WfTuTOU2b8EKTcPrlFzAlPYNeBYXuIwNnEZgb5+2QdRFmt6I1ZOcX4wiS1000HHJE0EA1sTXhus7X8X5uAMsPH+XOjEx+t1q4qUUzHoiMINlsgnkXw74VRpfqN1pHBHNB26YUntWxn/CNjz/+mIMHDxpdRoMkQdDQDZtJnLIyOSOLJYnHuC0jky9DghjTug3vmPMoev8P8MF1kLrf6ErrveFdY/j4zxcSFRpY/cwN0GOPPUa3bt2Ij4+nZ8+e/Pzz2X1IljVkyBDq+szP2X38v//++xw6dIg2bdrUaX2+cuutt7JjR/nOC8/VpEmTKu3iwxPk1FBDd8ZdR0GZiUx1hXF1tyk8nbWNZxO/ZVHHeO479hMX/vsC6D8ZBt0LNkfV6xT1n4fH0Pjxxx9ZtmwZGzduJDAwkJSUFAoLvdcB39l9/N94441e21ZliouLsVhq9xZ5do+l/kKOCBqD+HFw1zaYlQF3baNVv8m8NOwlXhr6EkW2UG6PDOXu83pw9JdX4F99YMM74HJWu9rG5reUk1zy7GpW7So/bkK9cmowpszDgD59p9g5jMx37NgxIiMjSztmi4yMpEULd3fcK1asoFevXvTo0YP/+7//o6CgoNzyZ47Q9emnnzJp0iQATpw4wdixY0lISCAhIaE0AE7Nr7Xm3nvvpXv37vTo0YOPP/4YcI9qNmTIEK655ho6d+7MDTfcUGEndK+99hr9+vUjISGBP/7xj+Tm5gLuT9iTJ09m4MCBdOzYkWXLlgHu8RCuvfZaxowZw2WXXUZaWhpXX3018fHx9O/fny1btlBcXEy/fv1YvXo1APfffz8PPvggUPYoKCQkhBkzZtCnTx8uvfRSfvnlF4YMGUK7du1YsmQJAAcPHmTgwIH07t2b3r17l/79WmumTJlC165dGTVqFElJp/e52bNn069fP7p3787tt99eZed7Naa19quvPn36aOE5+cX5ev7m+brve31133d763lvXKjzZ4Vp/crFWh9ca3R59Upieq5uPWOZ/uiX332+7R07dpz+ZfkMrd8cWfnX7CitHw4r/zU7qvJlls+ocvvZ2dk6ISFBd+jQQd9xxx169erVWmut8/LydGxsrN69e7fWWusJEyboOXPmaK21Hjx4sF63bp3WWuvg4ODSdS1cuFBPnDhRa631uHHjSucvLi7WGRkZZeb/9NNP9aWXXqqLi4v18ePHdVxcnD569KhetWqVDgsL04cPH9ZOp1P3799ff/fdd+XqTklJKf35wQcf1HPnztVaaz1x4kR9+eWXa6fTqffs2aNbtmyp8/Ly9FtvvaVbtmypU1NTtdZaT5kyRc+aNUtrrfWKFSt0QkKC1lrrbdu26c6dO+uvvvpK9+zZUxcUFJT7mwG9fPlyrbXWV199tR4+fLguLCzUmzZtKl3PyZMndV5entZa6z179uhT72+fffZZ6d995MgR7XA49MKFC7XWurQ2rbW+8cYb9ZIlS8r93WX2lxLAel3J+6ocETRygeZAbou/jSVXL2FQ3BBeMmdzdccEVhelod+6AhZOgoxDMu4z7ucIgPr/LIGz/CfyKqfXQEhICBs2bGD+/PlERUUxfvx43n77bXbv3k3btm1Lh6ucOHEia9asqfF6V65cyR133AHUbswBoHTMAZPJVDrmwNm2bdvGwIED6dGjBwsWLGD79u2lbePGjcNkMtGhQwfatWtX2lnd8OHDadq0aen2J0yYALh7YE1NTSUzM5Nu3boxYcIExowZw5tvvklAQEC5bQcEBDBixAgAevToweDBg7FarfTo0aO01qKiIm677TZ69OjBtddeW3p9Yc2aNaV/d4sWLRg69HT3M6tWreKCCy6gR48erFy5sszfVFdyjUAA0DykOc8NeY6fj/3MEz8/wV/D0hkY058Z+7+m9Y6l7qfYXCV9tDfSh9JCAi2YTYqsvNr3Ve9RV1Q9mDpzupecFjqLIw5u/m+dN2s2mxkyZAhDhgyhR48evPPOO6XDVVbnzPEHatNHvz7HMQcmTZrE559/TkJCAm+//Xbp6Zyzazrz9+Dg4Cq3f2q+rVu3Eh4ezokTJyqsz2q1ls5rMplK6zWZTKW1zpkzh5iYGDZv3ozL5cJmOz1mR0VjNuTn53PnnXeyfv164uLimDVrlkfGPJAjAlHGBc0vYOGVC7m3771sLM5kbMtoXmwaTq4+65pBI3woTSnlH08XD5vpHjPjTFa7e3od7d69m71795b+vmnTJlq3bk3nzp05ePAg+/btA+C9995j8ODB5ZaPiYlh586duFwuFi9efLrUkj7+AZxOJ1lZZQcxHDRoEB9//DFOp5Pk5GTWrFnD+eefX+O6s7Ozad68OUVFRSxYsKBM28KFC3G5XOzfv58DBw7QqVOncssPGjSodLnVq1cTGRlJWFgYixYtIjU1lTVr1jB16lQyMjJqXNOZMjMzad68OSaTiffeew+n01m63Y8++gin08mxY8dYtWoVcDpEIyMjycnJ8didRBIEohyrycpN3W5i2dhlXNF2JK+HBXFlbHO+DA4qO45oZvnB7hu6S7vEcF50JSPJ1Rfx42DMXPcRAMr9fczcczp6y8nJYeLEiXTt2pX4+Hh27NjBrFmzsNlsvPXWW1x77bX06NEDk8lUZoD7U5588klGjx7N0KFDad789JgdL774IqtWraJHjx706dOn3GmOsWPHEh8fT0JCAkOHDi0dc6CmHn30US644AKGDx9ebnyATp06MXjwYK644grmzZtX5tP4KbNmzWL9+vXEx8dz33338c4775CSksJ9993HG2+8QceOHZkyZQrTpk2rcU1nuvPOO3nnnXfo378/e/bsKT0aGTt2LB06dKBHjx7ccccdpeEaHh5eeirp6quvpl+/fnXa7tlkPAJRrU3/6s7jgUXsDAxgfFY296emYwb3G8xd5QcUEZ4n4xF41qRJkxg9ejTXXHON0aV4RW3HI5AjAlGtnoNn8mFyJjdnZPFxWCj3REeSrxQ0SwA/+yAhhChPgkBUL34c5jFzudsVyn2p6awIsvPn1u3J3Ptf+M8Ud59FjcRDn29lzL++N7oMcY7efvvtBns0UBcSBKJmSh5Ku+HuRJ4e8ixbzZqJ53Xn+NYP4cProSDH6Ap9wumCY5nnfpeGEPWJBIGotRFtRvDq8Fc5gZMb2nVi76Fv4Z3RkJNsdGle57Bbycor8szTnELUExIEok76NevH2yPeBmsQE1u1YV3mfnhjeIPvvC7MbqHQ6aKgWHogFQ2HBIGos05NO/H+yPeJCmnBn2Mi+YqT8MZlkLjB6NK8xmF3j0NQ758laICkG2rvkSAQ56R5SHPeveJdukV2555wOx+EBrlPE+35n9GleUWnmFCu7ROLqYKnPhs66Ya6rIMHD/LBBx8YXYZHSBcT4pw5Ah28dtlrzFgzgycOryQpIJZpH16PGvMC9L7J6PI8qm+bpvRt09ToMnyuMXZDXZ1TQfCnP/3J6FLOmRwRCI+wWWw8P+R5xnUcxxuWfB5q05miJX+F1U81uGcNtNaNbtxif+2Gev/+/YwYMYI+ffowcOBAdu3aVWU30pV18bxv3z4uvfRSEhIS6N27N/v37+e+++7ju+++o2fPnsyZM6fSLqX9gRwRCI8xm8w81P8hooOieWnTS6SeF8/z3z5BUNYRaDUAVv3TYwOlGOX31JMMe+5bnr4mnj/0jjWkhqd+eYpdabs8us7OTTsz4/wZlbZfdtllzJ49m44dO3LppZcyfvx4Bg8eTH5+PpMmTWLFihV07NiRm266iVdeeYXp06fXaLtTp05l8ODBLF68GKfTSU5O2duQFy1axKZNm9i8eTMpKSn069ePQYMGAfDrr7+yfft2WrRowUUXXcTatWu5+OKLyyx/++23M2/ePDp06MDPP//MnXfeycqVK0ufI5g7dy5ffvll6WmuKVOmMHOmu0+mCRMmsGzZMsaMGcMNN9zAfffdx9ixY8nPz8flcvHkk0/y7LPPlo5lkJuby9dff43NZmPv3r1cf/31dT415msSBMKjlFL8OeHPRAVFMfvH2dzcIZ5/b3qPyF/fA11yp40f914aEmih2KXJamQXi091Q/3dd9+xatUqxo8fz5NPPkmvXr3KdUP973//u8ZBsHLlSt59912gdt1Qh4WFlXZDDZR2Q31mEOTk5PDDDz9w7bXXlk47dbRyZjfSP/74Y2k30qtWreLpp58mNzeXtLQ0unXrxpAhQzhy5Ahjx44FqLBPInB3KT1lyhQ2bdqE2Wxmz549NXoN6gMJAuEVf+jwByLtkfxt9d+4qWUL5h07Tqszb7k81XupnwVBWOldQ8Z1RV3VJ3dv8rduqF0uF+Hh4WzatKnC5c/uRrqyLp5r+sxIVV1K13dyjUB4zaDYQbxx+RtkK82EFjFsP3vwDj/svdRqNhEUYCYrv3EdEfhjN9RhYWG0bduWhQsXAu5Q2bx5M0CF3UhX1sVzWFgYsbGxfP7554D7qCI3N5fQ0FCys7NLt1dZl9L+QIJAeFV8VDzvZZuwuzQ3N49mU+AZYeAw5hz7uXLYrY3uOQJ/7YZ6wYIFvPHGGyQkJNCtWzf+85//VNqNdFVdPL/33nvMnTuX+Ph4BgwYwPHjx4mPj8disZCQkMCcOXMq7VLaH0g31ML7tnxCyn+nMSHKQb4y8dHR48Rghqv+7XenhgD+tWIvcU2DuLpXS59tU7qhFrUh3VCL+id+HJGjXmRuroWTJsVdMVEUaA3RXY2urE7+OqyDT0NACG+TIBC+ET+ODlO38cTQuWwNDGB2TAz6o+shN83oymqt2OlqdKeGRMMmQSB8aljrYdyRcAdLbCYWkAWf/h84DR4MvpZmfLaVkS9+Z3QZQniMBIHwuckJkxkaN5Rnmzj46egPsOIRo0uqlTC7xZDnCPztep4wRl32EwkC4XMmZeLxgY/TNrw997SI5fAv/4atnxpdVo057FayC4px+rCbCZvNRmpqqoSBqJLWmtTU1Fo/wyAPlAlDBFuDmXvJXK7773VMjW3NgiVTCIrsCM3jjS6tWmE290NlWXlFNAkOqGZuz4iNjSUxMZHk5IY/+I84NzabrfSJ65qSIBCGiQuL45nBz3DH13fwUFQkz350A6bbV0NwhNGlVenUmARZ+b4LAqvVStu2bX2yLdH4yKkhYagBLQZwd9+7+TpQMd+UA59OqvcXj3vEOrh7eEeCA+VzlGgYJAiE4W7qehNj2o3h3+EhrDyxDr6eaXRJVeoYE8rUYR2IDAmsfmYh/IDXgkAp9aZSKkkpta2SdodSaqlSarNSartS6mZv1SLqN6UUMy+cSfeI7tzfrBn7N7wKmz8yuqxKOV2aoxl5ja6/IdFwefOI4G1gRBXtfwF2aK0TgCHAc0op35xwFfWOzWJjziVzsAeGM7VlKzKXTYejvxpdVoWSswsY8ORKlm0+ZnQpQniE14JAa70GqOqxUQ2EKnf/tCEl89bvk8PCq5oFN+OFS17gqAn+HhOJ86MbIaf+3SUTZndfG5AjAtFQGHmN4CWgC3AU2ApM01q7ql5ENHQ9o3vyUP+H+MGqeNGaBwsngbN+veHarWasZiXdTIgGw8gguBzYBLQAegIvKaXCKppRKXW7Umq9Umq93Efd8P2x4x+5rtN1vBUWzLKUjfC/B40uqQylFGG2xtcVtWi4jLz/7WbgSe1+VHKfUuo3oDPwy9kzaq3nA/PB3Q21T6sUhvj7+X9nX8Y+ZrGRtpvepFtxAexfUW/GPHbYrY1uuErRcBl5RHAIGAaglIoBOgEHDKxH1CNWk5XnhjxHRHAM05q3IGXTu+6xjtGnxzze8olh9U0d1oFr+vjnwDpCnM2bt49+CPwIdFJKJSqlblFKTVZKnRq+6FFggFJqK7ACmKG1TvFWPcL/NLU15cWhc8lULv4WE0mZz9+nxjw2yNW9WjKkU7Rh2xfCk7x2akhrfX017UeBy7y1fdEwdG7amUeTU7k3OpInIpowMzX9dKOBYx4nZeWTerKQLs0rvKwlhF+RJ4tFvTfCEsGtGZksDAtlWXDQ6QYDxzx+ccVebnj9Z8O2L4QnSRCI+m/YTKbkFJKQX8BTEU1IM5nAYndfMDZIWMkA9tIttGgIJAhE/Rc/DvOYucwqCCDHZOLJiCbQ5mLD7xpyujS5hU7DahDCUyQIhH+IH8d5U7dxe8+/8EVIMGuOroUjGw0r51RX1PIsgWgIJAiEX7m1x62c52jH7MgIcpb+1bAuq0sHp5FuJkQDIEEg/IrVbOWRix4lyWzihaIj8NPLhtTRq1U4L4zvSbOw2g0JKER9JEEg/E58VDw3dLmBj8NC2bj2aUg/6PMaWoTbubpXS8KDpMNc4f8kCIRf+muvv9IyKIaHm4ZRsHQ6+PjuncJiFz8fSOVoRp5PtyuEN0gQCL8UZA1i5oDZHLSaeTV9I2z91Kfbzy0sZvz8n/hy23GfblcIb5AgEH5rQMsBXNluDG85HOz+5n7IrWr4C88KtcldQ6LhkCAQfu3v588gLDCMh0PMFP/vIZ9t12xShAZa5K4h0SBIEAi/5gh0cP+F/2B7YAALflsCB7712bZPPV0shL+TIBB+7/LWlzOk5SBeatqEw8unuXsm9YEwu5WsPBldVfg/CQLh95RSPHThTCwWG49YTqK/fcYn2501pivTL+3gk20J4U0SBKJBiAmO4a5+9/Kz3cbiLa/BiR1e3+YF7SLo3tLh9e0I4W0SBKLBuKbjNfSJTODZJg6Sl/4FXN7tEG738Wy+3nHCq9sQwhckCESDYVImZl38TwrMFh4vPATr3/Tq9hauP8y0j3716jaE8AUJAtGgtHG04Y6eU/gmOIhvvn8cMo94bVsOu5XcQidFTpfXtiGEL0gQiAZnYveJdAlrx2PhdjKX3+217YSVdEWdJbeQCj8nQSAaHKvJyiODniTdbOH51F9g51KvbEfGJBANhQSBaJC6RHRhYreJLAoN4eev/w75mR7fRpjdAkgQCP+n/G3M1b59++r169cbXYbwA/nF+VyzeDTOzCMssp6HPe03yEx0D3o/bOY5D3WZmVvEvuRsOjcLIzjQ4qGqhfAOpdQGrXXfitrkiEA0WDaLjYcHPkGi1cLLGVsg8zCg3d+XToUtn5zT+h1BVvq0biohIPyeBIFo0Po168c1+S7edYSyPeCMQWSK8mDF7HNad0Gxk4XrD7PzWNY5VimEsSQIRIN39/GjRDqdzIxsSpmz+ZmJ57RereHeT7ewclfSOa1HCKNJEIgGLzSsJQ+mprMnMIAPwkJPNzhiz2m9NquZAItJbh8Vfk+CQDR8w2YytAguys3j1XAHGSYTWO3uC8bnKMxmlTEJhN+TIBANX/w4GDOXe4psnDQp5oU7YOi53zUE4LBb5PZR4fckCETjED+O86Zu45q2o/k4LITfEtd6ZLUOGZxGNAASBKJRufP8ewk0WXk++Sc4svGc1/fC+F48d23Pcy9MCANJEIhGJcIewW09bmV1cBA/f32v+9afc9AqIohmDpuHqhPCGBIEotG5Mf5WWlhCebYwEeeeL89pXesOpvH6dwc8VJkQxpAgEI1OoDmQ6f0fYFdgAEtXPwTOuo87/O3uZB5fvhN/66pFiDNJEIhGaUS7UcQHxzHXkkvuxrfqvB6H3YpLQ06BDGIv/JcEgWiUlFLcO+hxki0W3l43BwpP1mk90gOpaAgkCESj1TO6JyOi+/GWXXHiu2fqtA5H6eA0ckQg/JcEgWjUpg98FJfJzNzdH0BO7fsMCpPBaUQDIEEgGrWWIS25sf1YlgQHsv2b+2u9fO9WTVj34KX0a9PEC9UJ4RsSBKLRu/X8e2iqAnj2xHfopN21WtZmNRMVGojFLP+VhP+SvVc0eqEBodzZ807W2wJZ9fXfarVskdPF81/v4Yf9KV6qTgjv81oQKKXeVEolKaW2VTHPEKXUJqXUdqXUt96qRYjq/LH7RNpZHTyft4+i39bUeDmzUvxr5V5+2p/qxeqE8C5vHhG8DYyorFEpFQ68DFypte4GXOvFWoSoksVk4Z4Bs/jdauXjlX+vcdcTJpMizCYdzwn/5rUg0FqvAdKqmOVPwCKt9aGS+WWYJ2Goi1sPY0BIG15RWWRuXlDj5cKkK2rh54y8RtARaKKUWq2U2qCUusnAWoRAKcXfhjxDjsnMvJ+fguKCGi3nsFvJypfnCIT/MjIILEAfYBRwOfAPpVTHimZUSt2ulFqvlFqfnJzsyxpFI9MxojN/aHYhHwVqfv/h+Rot47BbZbhK4deMDIJE4Eut9UmtdQqwBkioaEat9XytdV+tdd+oqCifFikan78MepwAZWLOzncgL6Pa+d+Y2I9P/nyh9wsTwkuMDIL/AAOVUhalVBBwAbDTwHqEACAyKJJbz7uWFTYr61Y8UO38NqsZk0n5oDIhvMObt49+CPwIdFJKJSqlblFKTVZKTQbQWu8EvgS2AL8Ar2utK73VVAhfmnDBvTRTATxzbCWu9N+rnHflrhM8uHirjyoTwvO8edfQ9Vrr5lprq9Y6Vmv9htZ6ntZ63hnzPKO17qq17q61fsFbtQhRWzaLjem9p7EzwMqyr6ZXOe+Oo1ks+PkQ+UVO3xQnhIfJk8VCVOKKbjfSw9qEF3N2kXv450rnK+2BNF8uGAv/JEEgRCVMysS9Ax8jyWLhnZWVj28cVtoVtQSB8E8SBEJUoVfcQC4Lbc9brlSSnusAs8JhTnfY8knpPKe7opZnCYR/kiAQohrTWw6nWClesjkBDZmHYenU0jBw2K3YrWa5RiD8lgSBENWI++EVbszM5vOQYHYFuD/9U5QHK2YD0CsunJ2PjuCi8yINrFKIupMgEKI6mYncmplJmMvFC03Cy0wHd9cUQvgzCQIhquOIJcyluS0ji7VBdn6xBZZOB3C5NHd9vIllW44aWKQQdSdBIER1hs0Eq53rsrNpVlzMnKbhaIvdPR13V9RfbDvG5sMZxtYpRB1JEAhRnfhxMGYugWFxTEnPZFtgIF+3v8A9vYS74zm5a0j4JwkCIWoifhzctY3Rdx/mPAL5V+5eik6eHp7SYZfBaYT/kiAQohbMJjPTek3hoMXM4m/uKZ0uo5QJfyZBIEQtDe4xkd6mEOalriM3w90hXWwTO8GBZoMrE6JuJAiEqCWlFHf1f4hks4kF39wNwAvX9eL1if0MrkyIupEgEKIOenYYxSXWSN7M3kVGkvSeLvybBIEQdTRt4D/JVYrXVt7D0s1HueH1n3C6Ku6YToj6TIJAiDpqH3cRV9lb8WF+IieO/MjafankyCD2wg9VGgRKKYsvCxHCH9059FkUsP34HEDGJBD+qaojgl98VoUQfqpZVFf+FNaFFaYMOgRulFtIhV+qKgikJy0hauDWYc8TrCEyZrEEgfBLVZ3+iVJK3V1Zo9b6eS/UI4TfcTjiuD6sN6+ZfuXY75/AeXcaXZIQtVLVEYEZCAFCK/kSQpS47YoXiHZq/nP4DbTLZXQ5QtRKVUcEx7TWs31WiRB+zB7UlDtiL+WRYytYte5Fhl5wl9ElCVFjco1ACA/59MBVtCrWzN3xNsXFhUaXI0SNVRUEw3xWhRANwO6kIvqrIew3uVi69lGjyxGixioNAq11mi8LEcLfhdmtJAVOJN5p4t8HPie/IMfokoSoEXmyWAgPcditZBVopne7hRMm+HD1/UaXJESNSBAI4SGnBqfpd/5fudhp5fWjq8g8mWx0WUJUS4JACA/p3tJBh+gQUIrpfe4mW8GbK++pfkEhDCZBIISHzBjRmWeuTQCgU/wNjCKEBakbOJHxm8GVCVE1CQIhvEEppgx4GBfwysq/GV2NEFWSIBDCQz5Zf5iBT68kv8gJQMuOVzDeHMnirD0cOLHJ2OKEqIIEgRAeUuR0cTgtr0zHc7cNfhy71vxr9X0GViZE1SQIhPCQMJsVoEwQNG01gEm2VnyTf4TNL3aBWeEwpzts+cSgKoUoT4JACA9x2N1BkHVWV9Q3xQ0notjJnMBiNBoyD8PSqRIGot6QIBDCQ04FwdljEgSte5PJGZlssNv4zm5zTyzKgxXSp6OoHyQIhPCQqNBALu0SXRoIpTIT+WN2Dq2LiniuaROKz5guRH0gQSCEh7QIt/P6xH70bdO0bIMjFitwd1oGBwKsfBYaUjpdiPpAgkAIbxs2E6x2LsnNo19ePv9u4iA7IMg9XYh6QIJACA+65NnVPPnFrrIT48fBmLkoRxz3pKWTYTLxWscL3dOFqAckCITwoIIiJyk5BeUb4sfBXdvo+vejXFlk4v2cvSRmHvJ9gUJUQIJACA8KK+mBtFKWAP7a714s2sULq+72XWFCVMFrQaCUelMplaSU2lbNfP2UUk6l1DXeqkUIXwmzW8s9R3C2mIQbuVmH8b/M3Ww68qOPKhOict48IngbGFHVDEopM/AU8D8v1iGEzziqOyIAUIqJlzxFdHExT695AJd2+aY4ISrhtSDQWq8Bqhvu8q/AZ0CSt+oQwpcGdYhkaOfoaucLajuIv9rasrUwhS93yRPGwliGXSNQSrUExgLzjKpBCE+bcGEb/j6ic43mvfKyF+hSUMgL658jvzjfy5UJUTkjLxa/AMzQWjurm1EpdbtSar1San1ysgz9J+o3l0ujta52PlNUR+6JGsAxVz7vr3/RB5UJUTEjg6Av8JFS6iBwDfCyUurqimbUWs/XWvfVWveNioryYYlC1M57P/1O+weXk55bzXWCEucPf5pL8gp5ffcHpOSleLk6ISpmWBBordtqrdtordsAnwJ3aq0/N6oeITwhyGpG6/I9kFYqJIq7O4yjQDt5+buHvVucEJXw5u2jHwI/Ap2UUolKqVuUUpOVUpO9tU0hjFZZD6RVaTPwfsbnu/js6Br2pu3xVmlCVMqbdw1dr7VurrW2aq1jtdZvaK3naa3LXRzWWk/SWn/qrVqE8JWwU2MS5Nc8CAgIYnLv6QS7XDz37QwvVSZE5eTJYiE8qC5HBADhfW9hcrGdtVn7+P7Qai9UJkTlJAiE8KDo0EAmDWhDq6ZBtVvQZOb6wf+kVVERz66dSbGruPplhPAQCQIhPKhJcACzruxGfGx4rZe1drycu60t2V+YzqIdCzxfnBCVkCAQwsPyi5ycLKjDJ3qlGDr8Ofrk5fPvX/9FTmGO54sTogISBEJ42MVPreSx5TvrtKxq0ZN7I84nzVXA6+vneLgyISomQSCEh4XZqu+BtCrdhj/BmJN5vLd3IUdyjniwMiEqJkEghIdVOyZBdcJbMbXd1ZhcTl5c+4jnChOiEhIEQnhYTcYkqE6zwQ8xMbeIL47/yKakTZ4pTIhKSBAI4WEOu5Ws/HO8/dMezv/1/AuRxU6e+f6hGnViJ0RdSRAI4WGj45sz8cLW57yeoAvuYGqhhS3Zv/O/377wQGVCVEyCQAgPu7xbMyZd1PbcV2QJ4MqLZ9KpoJAXfnqcAmfBua9TiApIEAjhYflFTg6n5VLsPPchKM3d/8g9piiOFGWy4NXeMCsc5nSHLTKqmfAcCQIhPOzzX48w8OlVJGV74BO8UvRvN4LBuXnMt0OqSUHmYVg6VcJAeIwEgRAeVqceSKuy+SPuTkunQCmea9rEPa0oD1bM9sz6RaMnQSCEh5X2QFrDUcqqlZlIu6Ji/i8ji6WhwXxvt5VOF8ITJAiE8LAwW926oq6UIxaAP2dk0rawiNmRTTmpVOl0Ic6VBIEQHuYoPTXkoa6kh80Eq50A4JGUVI6bzcyNiHBPF8IDJAiE8LCo0EBmju5KQqzDMyuMHwdj5oIjjl4FRVyXlcOHoUFsatbRM+sXjZ7ytycW+/btq9evX290GUIY5uSi27g6/XuCHK1ZOHYJAeYAo0sSfkAptUFr3beiNjkiEMIL9ifncDgt1yvrDh7xFDOzizlw8givbXnVK9sQjYsEgRBecNMbvzDn6z3eWXlQUwYOfYxROSd5fcvr7E3f653tiEZDgkAILzjnrqir0+NaZoR2J9RZzKzvHsDpcnpvW6LBkyAQwgscdovnHiiriFI0Gf0iMzJy2JK+iw92feC9bYkGT4JACC9wePuIAKBpW0aefxcDc/P414Y5JGbLA2aibiQIhPAC93CVHnqOoArqwin8wxSNKi7kkbX/kHELRJ1IEAjhBded34pZV3bz/obMFpqP/hfT0zP46cR6/rP/P97fpmhwJAiE8II+rZswonsz32ysZR/Gd7mRXvkFPPPzE6Tkpfhmu6LBkCAQwgtScgpYuy+F/CLf3M1jGvoQs/IDyCvK5YmfHvPJNkXDIUEghBd8vzeFG17/maMZeb7ZYGAI7a54nsnpGXx16BtWHlrpm+2KBkGCQAgvCLNbAA/2QFoTHS/j5pZD6VhYxGM/PEJ2Ybbvti38mgSBEF7g8R5Ia8h6xdPMziogJT+N59c/79NtC/8lQSCEF5QOTuPLIwKA0Bi6DZnFhKwsPt37KeuOr/Pt9oVfkiAQwgtODU6T5esgAOg1gb+EdiG22MWs7/9BfnG+72sQfkWCQAgvaBIcwGs39eWSztG+37jJhH30XB5OzeDQySO8vPll39cg/IoEgRBeYDWbGN41hpbhdmMKiOpI//OnMjY7h3e3vcOO1B3G1CH8ggSBEF6ydl8KGw+lG1fAxdP5G5E0cbl4+Pt/UOQy4DSV8AsSBEJ4ycNLtvP6dweMK8ASiGPMXB5ITmZXxh7e2f6OcbWIek2CQAgv8UkPpNVpfSHDu1zHsJO5vLLhBQ48FglzusOWT4ytS9QrEgRCeEmYzeKTHkir1bwnD6amEeRycVd0JDlZibB0qoSBKCVBIISX1IsjAoDvniXK6eK5pBR+t1q4PyoCV1EerJhtdGWinvBaECil3lRKJSmltlXSfoNSakvJ1w9KqQRv1SKEEbw+XGVNZboHrDk/v4B709JZHRzEy+GO0ulCePOI4G1gRBXtvwGDtdbxwKPAfC/WIoTP3XJxW9675XyjywBHbOmPf8rKYWx2Dq82cfB1ZGwVC4nGxGtBoLVeA6RV0f6D1vrUvXU/AbJXigaldUQw8bHhRpcBw2aC1f08gwIeSkkjPr+ABx0B7E7bbWxtol6oL9cIbgG+MLoIITzpcFouH/1yyPjTQ/HjYMxccMQBioDgaF5ISiHUpZm2ahoZ+RnG1icMZ3gQKKUuwR0EM6qY53al1Hql1Prk5GTfFSfEOdh+NJP7Fm0lMT3X6FLcYXDXNpiVAffuJWrE08w5mkhSzjHu+fYeil314O4mYRhDg0ApFQ+8DlyltU6tbD6t9XytdV+tdd+oqCjfFSjEOQgzqgfSmuh7C/Fdx/NwcjI/H/+Z59Y/Z3RFwkCGBYFSqhWwCJigtd5jVB1CeMvpHkjr4adtpWDUc1wV3pUbs/N4f+f7fL7vc6OrEgbx5u2jHwI/Ap2UUolKqVuUUpOVUpNLZpkJRAAvK6U2KaXWe6sWIYxQOjhNfTwiALAEwrj3+Fu+iQuKYPaPs9mSvMXoqoQBvHnX0PVa6+Zaa6vWOlZr/YbWep7Wel5J+61a6yZa654lX329VYsQRnAE1eNTQ6eENccy7j2ePX6caJdm+qrpJOUmGV2V8DHDLxYL0VCFBFj46q5BjOsXZ3QpVWt1AeFXPMPcxMPk5Kdx1+q7KHQWGl2V8CEJAiG8xGRSdIwJLT1FVK/1mUTH+Bt57PhxtiRv4dGfHkVrbXRVwkckCITwokUbE/ly23Gjy6iZK55meEQCf87K5fN9n/PBrg+Mrkj4iASBEF705trf+GT9YaPLqBlLAIx7lzsLrQwphGfWPcMvx34xuirhAxIEQnhRvemBtKZCYzCNf58nTpygtTbxt2//RmK2dE7X0EkQCOFFYTY/CwKA2L6EjHyOuYd/x1l4kmmrppFbVA+ejhZeI0EghBc57Nb6+xxBVXpPoHWvm3nm2BH2pe/lH2v/IRePGzAJAiG8yO9ODZ1pxBNcFN2HuzKy+er3r3h96+tGVyS8RPlbyvft21evXy8PIQv/kJlXhMulaRIcYHQpdZOTjJ4/mPuDXCy3WXg8OY3Rlqburq3jxxldnagFpdSGyh7clSMCIbzIYbf6bwgAhESh+kxi1onj9Mkv4P7oCN7RGTLmcQMjQSCEF+1LyuHpL3eRlJVvdCl1t/FdbFoz70QSw0/m8mxEE54JDcQlYx43GBIEQnhRYnouL6/ez+H0PKNLqbuSsY0DNTyTlMKfMrN51xHGfQG50hVFAyFBIIQX1fseSGvijDGPzcB9aelMT0vni5Bg7vzmTnIKc4yrTXiEBIEQXlSvB6epqTPGPAb3uMe3ZGbzeHoeG06sZ9KXk0jOlZED/ZkEgRBeVHpEkO/HQXDWmMc44mDYw4xxWngpJZNDmQe5cfmN/Jb5m9GVijqSIBDCi06PUubHQQBlxzy+axsMvBtu+ZqLAqN568gR8gsymfDFBDYlbTK6UlEHEgRCeFGAxcTO2SP4yyXnGV2K54XHwf99Sbfonrz/214cLs1tX93G6sOrja5M1JIEgRBeZg8wo5QyugzvsDeBCYuJ6zCSd/ft4DyTjWmrpvHpnk+NrkzUggSBEF4279v9vPPDQaPL8B6rDa59m4i+t/HG3q0MMIXyyI+P8PKml6V/Ij8hQSCEl32z4wRfbDtmdBneZTLDFU8RdOkjzN23lat0MK9sfoVHfnyEYlex0dWJakgQCOFl7h5IG8GboVJw0TSsY+fz6KG93FYYwGd7P+OuVXeRV+zHD9Q1AhIEQniZX/dAWhcJ41E3LGRq8nEezHHxbeK33PrVraTnpxtdmaiEBIEQXhZmt/r3cwR10X4o3Lyc6/KKeD49h13J27jpg4EceCwS5nSXDuvqGQkCIbzMYbeiAJerkV04bZ4At3zNpUUW5h89Rjou/tiyGc+Zsjm5bJqEQT0i4xEI4WVa64Z7+2hNPN8Vso6QajIxt2k4i0JDiCou5u58M6Pu3Nq4XxsfkvEIhDBQo3+jyzoKQITLxSMpaXxw5DgxTif3hygmfTmJXWm7DC5QSBAI4WU7jmYx5YON/JZy0uhSjHFG76UAPQoLWXD0BI+kZvBb+h7GLxvPP3/6J5kFmQYVKCQIhPCyrPwilm05xtGMRnoL5Vm9lwKYzAH8odDE0n07GW+OZOHuhYxePJpP93yK0+U0qNDGS4JACC9rEGMSnIuKei+96t9w904cA2fwwG/b+eR4Cu2w8siPj/Cn5X9ic/Jmo6tuVCxGFyBEQ+doCGMSnKv4cRUPdn/J/dDrRjp98zBvb/uM5ZEted56mBuX38jV513NtN7TiLRH+r7eRkaOCITwsgYxOI03hcfBNW+ibv6SUdYoluzdwc2uYJbtX8qYxWN4f8f70k2Fl0kQCOFlwQFmYsICMZsa+d1D1Wl9Idy2iuArX+LulFQ+O3yYeG3hqXVPce3Cy1j3r+4wK1weSPMCeY5ACFH/5GfBd8+hf3qZlXYbTzcJ4ajFQq/8fEbl5HJ5oYvwUS9WfLpJVKiq5wgkCIQQ9VfqfnjlIvKd+XwYFsJ/QoLZHxCARWsuLoJRw59lSOwQbBab0ZXWe1UFgVwsFsIHnli+E4D7R3YxuBI/E9EeivOxobk5M5tJmdnsDrCyLCSY5cFBrP72XoKtwQxvPZxR7UbRL6YfZpPZ6Kr9jgSBED6w41gW2flywbNOHLGQeRgABXQuLKJzWgZ3pWWwzmZjWbO2fH3gv3y+73Oi7dGMbDeSUe1G0alJJ3mqu4YkCITwgTC7lSON9YGyczVsJiydCkVnvH5WO+bB99G/OJ/+2xbxUOoBVgcF8d9oO+9vf5e3t7/NeeHnMardKEa2HUmLA9/DitmQmegOlmEz5frCGSQIhPAB9+A0cvtonZx6w67sjXzwDGxJOxmxfTEjti8iPf0QX4WEsIzDvLjxRV7c+CJ9CgoZ5cqmv8VEbOZh1NKpZdfdyEkQCOEDYTb34DSNvifSuqrsgTRwj4wW09X9dckDNDmxnfHbFzN++yIOZx1heUgwy4KDmR0ZAUCo00XXwkK6/jCbrqFhdInoQlxoHCbVeO+mlyAQwgfimtrpEB1KodNFoEUuZnqNUtCsu/tr6EPEHd/Cn18dxO0ZWewJsLIlMIAdAQHsCAzg/QAoWnMvACGWYDpHdKFrRFe6lHxvHdr69IXnLZ806FNLcvuoEKJhm9O99GLzmYqUiX1WMztLgmGnPZjdVjMFuN8T7eZAujTtQhcC6Lr7G7rkniS2uBi71u5O9MbM9aswMOQ5AqXUm8BoIElr3b2CdgW8CIwEcoFJWuuN1a1XgkAIUStbPqnwYjNj5kLbwZC0HU7sgBPbKUraxoH0/ey0wI6AAHYGBrA7wEqe6fRpo1CniyinkyjMRHYcRXRwMyLtkUQHRZf5HmQNKlvDuRxReOCIxKggGATkAO9WEgQjgb/iDoILgBe11hdUt14JAuGPth3J5MHFW/nn1T3oEeswupzGpzZvpM5iSDtQEhDbca55hoNWCzsDAjhusZBsNpNsMbu/l/xcWMF1n2CTlagAB1FaEZV2mKjiQkJdLoJcmiBlwd7rJoLaDyXIEoTdYifIWvLdEkSQNQiryeq+nlRVkNUiDAx5oExrvUYp1aaKWa7CHRIa+EkpFa6Uaq61PuatmoQwitOl2ZyYydc7jnM8Kx+LWXFJp2gANh/OICm7oMz8NquJgR2iANjwezppJwvLtAcHmhnQ3t0r5y+/pZXr0M5ht3J+26YA/LA/hZMFZfv4bxocQJ/WTQD4fm8KeUVl26NDA0mICwdg9e4kipxlPzA2d9jo3tIdaCt2nuDs4Zhjm9jp0jwMl0uzYldSudejTUQQHWJCKSx28e2e5HLt7aOCaRcVQl6hk+/3pZRr7xQTSquIILLzi/jpQFq59q4twmgZbicjt5B1B9PBMhguX1HaHh/rIAZIySng10MZ5Zbv1ao1kd06ciLuCsLWL6B97lHaF531HIi9CVk9byc96QgFhUlkF6aSWZxBhjObdJ1fEhQZJFvMbAkwk2wPoeCMIwt+X+z+qoQZhU1ZCC4qICg6HLt2MDY7h+uzc9yhsGK2x05NGXmxuCVw5om7xJJp5YJAKXU7cDtAq1atfFKcEJ7UzGFDKZi7ch8A4UFWNs28DIB53+7ni23Hy8zfMtzO2vuGAvDCN3v4bm/ZN8OOMSF8dddgAJ78Yicbz3oz69UqnMV3XgTAI0t2sPtEdpn2gR0iee8W9wH4fYu2kJhe9hmHEd2aMW9CHwDu+ngT6bllg+YPvVvy/LieANzx/kYKna4y7Tdd2JrZV3Wn2KW57d3yR/B3DGnPjBGdOVlQXGH7vZd34i+XnEfqyYIK22eN6cqki9pyNCO/wvanr4lnXN849iefrLD95Rt6M7JHc3Yczaqw/Z3/O5/BHaPY+Hs6X2RczZPW1wlSp8PYabZjvuJpvioawD2ryo+d8OXUAXQOK+I/azcz5odrONXfYBGQa1LkKRO5JsXeJhdwOD0dzAUoVYTLVITLVIzJ6iJPF5OnFLkmE7lKkWcyYTvzDE5mYrnt1pVXLxaXHBEsq+TU0H+BJ7TW35f8vgL4u9Z6Q1XrlFNDwl8dTsst/eRuUoquLcIAOJSaS1Z+2Tdaq9lEp2ahAPyWcpKTBWU/jQZaTHSIcbfvT84hr7DsJ3p7gJn2USEA7D2RTUFx2Tfq4EALbSODAdh9PJuis97Iw2xWWkW4z3HvOJqF66z3CYfdSlxTd/u2I+WHmGwaHECLcDsul2bHsaxy7ZEhgTRz2Ch2uth1PLtce3RYINGhNgqKnew9kVOuvZnDRmRIIPlFTvYllW9vGW6nSXAAJwuKKxwiNK5JEI4gK9n5RfyemluuvXVEEKE2K5m5RRxOz8WxbzEx657GmnOUopAWuIbOxNb7OtJPFlb4oOB50SHYrGZScgoIm9eLgJwj5ebBEUfSretJyioo19S5WSgWBSfS0mny5gACco9XuDx3bSs/vRKGdTpXTRC8CqzWWn9Y8vtuYEh1p4YkCIQQfuVcz/H74BqBkU9QLAFuUm79gUy5PiCEaHAqGqqzNm/i57p8DXjzrqEPgSFAJHACeBiwAmit55XcPvoSMAL37aM3a62r/agvRwRCCFF7Rt01dH017Rr4i7e2L4QQomYab+caQgghAAkCIYRo9CQIhBCikZMgEEKIRk6CQAghGjkJAiGEaOQkCIQQopHzu4FplFLJwO91XDwSKN+VofHqa11Qf2uTumpH6qqdhlhXa611VEUNfhcE50Iptb6yJ+uMVF/rgvpbm9RVO1JX7TS2uuTUkBBCNHISBEII0cg1tiCYb3QBlaivdUH9rU3qqh2pq3YaVV2N6hqBEEKI8hrbEYEQQoizNJggUEqNUErtVkrtU0rdV0G7UkrNLWnfopTqXdNlvVzXDSX1bFFK/aCUSjij7aBSaqtSapNSyqODMNSgriFKqcySbW9SSs2s6bJeruveM2rappRyKqWalrR58/V6UymVpJSqcGxAA/ev6uoyav+qri6j9q/q6vL5/qWUilNKrVJK7VRKbVdKTatgHu/uX1prv/8CzMB+oB0QAGwGup41z0jgC0AB/YGfa7qsl+saADQp+fmKU3WV/H4QiDTo9RqCe5jRWi/rzbrOmn8MsNLbr1fJugcBvYFtlbT7fP+qYV0+379qWJfP96+a1GXE/gU0B3qX/BwK7PH1+1dDOSI4H9intT6gtS4EPgKuOmueq4B3tdtPQLhSqnkNl/VaXVrrH7TW6SW//gTEemjb51SXl5b19LqvBz700LarpLVeA6RVMYsR+1e1dRm0f9Xk9aqMoa/XWXyyf2mtj2mtN5b8nA3sBFqeNZtX96+GEgQtgcNn/J5I+Reysnlqsqw36zrTLbhT/xQNfKWU2qCUut1DNdWmrguVUpuVUl8opbrVcllv1oVSKgj3MKefnTHZW69XTRixf9WWr/avmvL1/lVjRu1fSqk2QC/g57OavLp/eW2oSh9TFUw7+3aoyuapybJ1VeN1K6Uuwf0f9eIzJl+ktT6qlIoGvlZK7Sr5ROOLujbifiQ9Ryk1Evgc6FDDZb1Z1yljgLVa6zM/3Xnr9aoJI/avGvPx/lUTRuxfteHz/UspFYI7eKZrrbPObq5gEY/tXw3liCARiDvj91jgaA3nqcmy3qwLpVQ88DpwldY69dR0rfXRku9JwGLch4E+qUtrnaW1zin5eTlgVUpF1mRZb9Z1hus467Ddi69XTRixf9WIAftXtQzav2rDp/uXUsqKOwQWaK0XVTCLd/cvT1/4MOIL95HNAaAtpy+YdDtrnlGUvdjyS02X9XJdrYB9wICzpgcDoWf8/AMwwod1NeP0cybnA4dKXjtDX6+S+Ry4z/MG++L1OmMbbaj84qfP968a1uXz/auGdfl8/6pJXUbsXyV/97vAC1XM49X9q0GcGtJaFyulpgD/w30V/U2t9Xal1OSS9nnActxX3vcBucDNVS3rw7pmAhHAy0opgGLt7lQqBlhcMs0CfKC1/tKHdV0D3KGUKgbygOu0e88z+vUCGAt8pbU+ecbiXnu9AJRSH+K+0yVSKZUIPAxYz6jL5/tXDevy+f5Vw7p8vn/VsC7w/f51ETAB2KqU2lQy7QHcIe6T/UueLBZCiEauoVwjEEIIUUcSBEII0chJEAghRCMnQSCEEI2cBIEQQjRyEgRCCNHISRAIIUQjJ0EgxDlSSvUr6SPeppQKLulTvrvRdQlRU/JAmRAeoJT6J2AD7ECi1voJg0sSosYkCITwAKVUALAOyMfdr4/T4JKEqDE5NSSEZzQFQnCPMGUzuBYhakWOCITwAKXUEtyjQ7UFmmutpxhckhA11iB6HxXCSEqpm3D36vmBUsoM/KCUGqq1Xml0bULUhBwRCCFEIyfXCIQQopGTIBBCiEZOgkAIIRo5CQIhhGjkJAiEEKKRkyAQQohGToJACCEaOQkCIYRo5P4fhw8EeUq8Qd4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "T_exacta = np.zeros(N)\n", "error = np.zeros(N)\n", "\n", "print(b)\n", "T_exacta = 1.5-0.5*special.erf((x-1)/(2*np.sqrt(b*tf)))\n", "\n", "fig, ax1 = plt.subplots(figsize=(6, 5))\n", "ax1.plot(x,T[0, :],'--',label='Condicion inicial')\n", "ax1.plot(x,T[nt-1,:],'-o',label=f'Solución aproximada')\n", "ax1.plot(x,T_exacta,'-',label='Solución exacta')\n", "ax1.set_xlabel(\"x\") \n", "ax1.set_ylabel(\"T\") \n", "ax1.legend()" ] } ], "metadata": { "colab": { "provenance": [], "toc_visible": true }, "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.6" } }, "nbformat": 4, "nbformat_minor": 1 }