{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Ecuaciones en derivadas parciales con Python" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Esta notebook fue creada originalmente como un blog post por [Raúl E. López Briega](https://relopezbriega.com.ar/) en [Matemáticas, análisis de datos y python](https://relopezbriega.github.io). El contenido esta bajo la licencia BSD.*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\"Ecuaciones" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introducción\n", "\n", "Nuestra comprensión de los procesos fundamentales de la naturaleza se basa en gran medida en [Ecuaciones en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales). Ejemplos de ello son las vibraciones de los sólidos, la dinámica de los fluidos, la difusión de los productos químicos, la propagación del calor, la estructura de las moléculas, las interacciones entre fotones y electrones, y la radiación de ondas electromagnéticas. Las [Ecuaciones en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) también juegan un papel central en las matemáticas modernas, especialmente en la [geometría](https://es.wikipedia.org/wiki/Geometr%C3%ADa) y el [análisis](https://es.wikipedia.org/wiki/An%C3%A1lisis_matem%C3%A1tico); lo que las convierte en una herramienta de suma utilidad que debemos conocer.\n", "\n", "**_Nota_**: Este artículo corresponde a la tercer entrega de mi serie de artículos sobre [Cálculo](https://es.wikipedia.org/wiki/C%C3%A1lculo_infinitesimal) con [Python](https://www.python.org/); los anteriores fueron: [Introducción al Cálculo](https://relopezbriega.github.io/blog/2015/12/02/introduccion-al-calculo-con-python/) y [Ecuaciones Diferenciales con Python](https://relopezbriega.github.io/blog/2016/01/10/ecuaciones-diferenciales-con-python/), los cuales es recomendable haber leído previamente.\n", "\n", "## ¿Qué es una ecuación en derivadas parciales?\n", "\n", "Una [Ecuación en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) es una ecuación que, como su nombre lo indica, contiene [derivadas parciales](https://es.wikipedia.org/wiki/Derivada_parcial). A diferencia de lo que habíamos visto con las [ecuaciones diferenciales ordinarias](https://relopezbriega.github.io/blog/2016/01/10/ecuaciones-diferenciales-con-python/), en donde la función incógnita depende **solo** de *una variable*; en las [Ecuaciones en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales), o [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) para abreviar, la función incógnita va a depender de **dos o más** *variables independientes* $x, y, \\dots$. Generalmente a la función incógnita la vamos a expresar como $u(x, y, \\dots)$ y a sus derivadas parciales como $\\partial u / \\partial x = u_x$ o $\\partial u / \\partial y = u_y$ dependiendo de sobre que variable estemos derivando. Entonces una [Ecuación en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) va a ser la identidad que relaciona a las *variables independientes* ($x, y, \\dots$), con la *variable dependiente* $u$ (nuestra función incógnita), y las [derivadas parciales](https://es.wikipedia.org/wiki/Derivada_parcial) de $u$. Lo podríamos expresar de la siguiente forma:\n", "\n", "$$F(x, y, u(x, y), u_x(x, y), u_y(x, y)) = F(x, y, u, u_x, u_y) = 0$$\n", "\n", "Al igual de como pasaba con las [ecuaciones diferenciales ordinarias](https://relopezbriega.github.io/blog/2016/01/10/ecuaciones-diferenciales-con-python/), el ***orden*** de una [Ecuación en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) va a estar dado por la mayor [derivada](https://es.wikipedia.org/wiki/Derivada) presente. Por lo tanto, el caso que expresamos más arriba corresponde a una [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales), con dos *variables independientes*, de *primer orden*. Si quisiéramos expresar una [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) de *segundo orden*, podríamos hacerlo de la siguiente manera:\n", "\n", "$$F(x, y, u, u_x, u_y, u_{xx}, u_{xy}, u_{yy})=0.$$\n", "\n", "### Clasificación de ecuaciones en derivadas parciales\n", "\n", "La clasificación de las [Ecuaciones en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) va a ser algo fundamental, ya que la teoría y los métodos para poder solucionarlas van a depender de la *clase* de ecuación con la que estemos tratando. Las clasificaciones más importantes que debemos tener en cuenta, son:\n", "\n", "1- **El *orden* de la [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales):** Como ya mencionamos, el mismo va a estar dado por el *orden* de la mayor [derivada](https://es.wikipedia.org/wiki/Derivada) presente.\n", "\n", "\n", "2- **Número de variables:** Esta clasificación va a estar dada por la cantidad de *variables independientes* que contenga la [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales).\n", "\n", "\n", "3- **Linealidad:** Esta es una de las clasificaciones más importantes, vamos a poder clasificar a las [Ecuaciones en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) en [lineales](https://es.wikipedia.org/wiki/Lineal) o [no lineales](https://es.wikipedia.org/wiki/No_linealidad). En las [lineales](https://es.wikipedia.org/wiki/Lineal), la *variable dependiente* $u$ y todas sus derivadas, van a aparecer en una forma *lineal*, es decir, que van a tener *grado* uno (no van a estar elevadas al cuadrado, o multiplicadas entre sí). Más precisamente, una [Ecuación en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) de segundo orden con dos variables, va a tomar la siguiente forma:\n", "\n", "$$Au_{xx} + Bu_{xy} + Cu_{yy} + Du_x + Eu_y + Fu = G$$\n", "\n", "en donde $A, B, C, D, E, F,$ y $G$ pueden ser *constantes* o una *función* dada de $x$ e $y$. Por ejemplo:\n", "\n", "$u_{tt} = e^tu_{xx} + \\sin t$, sería una ecuación [lineal](https://es.wikipedia.org/wiki/Lineal).\n", "\n", "$uu_{xx} + u_y = 0$, sería una ecuación [no lineal](https://es.wikipedia.org/wiki/No_linealidad).\n", "\n", "$u_{xx} + yu_{yy} + u_x = 0$, sería una ecuación [lineal](https://es.wikipedia.org/wiki/Lineal).\n", "\n", "$xu_x + yu_y + u^2 = 0$, sería una ecuación [no lineal](https://es.wikipedia.org/wiki/No_linealidad).\n", "\n", "**Tipos de ecuaciones lineales:** Asimismo, a las [Ecuaciones en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) [lineales](https://es.wikipedia.org/wiki/Lineal) de segundo orden las vamos a poder subdividir en las siguientes categorías:\n", "\n", "* **[Ecuaciones parabólicas](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_parab%C3%B3lica_en_derivadas_parciales)**: Las cuales van a describir el *[flujo del calor](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_del_calor)* y el proceso de *[difusión](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_de_difusi%C3%B3n)*. Éstas ecuaciones van a satisfacer la condición $B^2 -4AC = 0$.\n", "\n", "* **[Ecuaciones hiperbólicas](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_hiperb%C3%B3lica_en_derivadas_parciales)**: Las cuales describen los *[sistemas de vibración](https://es.wikipedia.org/wiki/Vibraci%C3%B3n)* y los *[movimientos de ondas](https://es.wikipedia.org/wiki/Onda)*. Satisfacen la condición $B^2 -4AC > 0$.\n", "\n", "* **[Ecuaciones Elípticas](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_el%C3%ADptica_en_derivadas_parciales)**: Las cuales describen los fenómenos de *[estados estacionarios](https://es.wikipedia.org/wiki/Estado_estacionario)* y satisfacen la condición $B^2 -4AC < 0$.\n", "\n", "4- **Homogeneidad:** Otra clasificación que podemos utilizar, es la de homogeneidad. Una ecuación va a ser *homogénea*, si el lado derecho de la ecuación, $G(x, y)$, es idénticamente cero para todo $x$ e $y$. En caso contrario, se la llama *no homogénea*.\n", "\n", "5- **Tipos de coeficientes:** Por último, podemos clasificar a las [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) de acuerdo a sus coeficientes $A, B, C, D, E, $ y $F$, si los mismos son *constantes*, se dice que la ecuación es de *coeficientes constantes*, en caso contrario será de *coeficientes variables*.\n", "\n", "### ¿Cómo resolver ecuaciones en derivadas parciales?\n", "\n", "Existen varios métodos que podemos utilizar para intentar resolver las [Ecuaciones en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales), la principal idea detrás la mayoría de estos métodos es la transformar a estas ecuaciones en [ecuaciones diferenciales ordinarias](https://relopezbriega.github.io/blog/2016/01/10/ecuaciones-diferenciales-con-python/) ([EDO](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria)), o en alguna [ecuación algebraica](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_algebraica); las cuales son más sencillas de resolver. Algunos los métodos que podemos utilizar son:\n", "\n", "1- **El [método de separación de variables](https://es.wikipedia.org/wiki/M%C3%A9todo_de_separaci%C3%B3n_de_variables):** Este método es uno de los más importantes y más productivos a la hora de resolver a las [Ecuaciones en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales). La idea es reducir a la [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) de $n$ variables, en $n$ [EDOs](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria).\n", "\n", "2- **El método de la [transformada integral](https://es.wikipedia.org/wiki/Transformada_integral):** Este método es similar al que ya vimos al resolver [EDOs](https://relopezbriega.github.io/blog/2016/01/10/ecuaciones-diferenciales-con-python/). La idea es aplicar una [transformada integral](https://es.wikipedia.org/wiki/Transformada_integral) para reducir una [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) de $n$ variables, en otra de $n - 1$ variables. De esta forma una [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) de 2 variables, puede ser transformada en una [EDO](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria).\n", "\n", "3- **El método del [cambio de coordenadas](https://es.wikipedia.org/wiki/Sistema_de_coordenadas#Cambios_de_coordenadas):** Este método intenta cambiar a la [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) original en una [EDO](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) o en otra [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) más sencilla de resolver por medio del cambio de coordenadas del problema.\n", "\n", "4- **El método de [perturbación](https://es.wikipedia.org/wiki/Teor%C3%ADa_perturbacional):** Este método aplica la [teoría perturbacional](https://es.wikipedia.org/wiki/Teor%C3%ADa_perturbacional) para intentar cambiar un problema de [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) [no lineal](https://es.wikipedia.org/wiki/No_linealidad) en una serie de problemas [lineales](https://es.wikipedia.org/wiki/Lineal) que se aproximan al [no lineal](https://es.wikipedia.org/wiki/No_linealidad) original.\n", "\n", "5- **El método de expansión de [autofunciones](https://es.wikipedia.org/wiki/Autofunci%C3%B3n):** Este método intentan encontrar la solución de una [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) como una suma infinita de [autofunciones](https://es.wikipedia.org/wiki/Autofunci%C3%B3n). Estas [autofunciones](https://es.wikipedia.org/wiki/Autofunci%C3%B3n) son halladas por medio de la resolución del problema del [valor propio](https://es.wikipedia.org/wiki/Vector_propio_y_valor_propio) que corresponde al problema original.\n", "\n", "6- **El método de las [ecuaciones integrales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_integral)**: Este método convierte a la [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) en una [ecuación integral](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_integral); la cual luego es resuelta aplicando ciertas técnicas particulares que se aplican a ese tipo ecuaciones.\n", "\n", "7- **Métodos numéricos**: La mayoría de las técnicas para resolver numéricamente a las [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) se basan en la idea\n", "de discretizar el problema en cada variable independiente que se produce en la [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales), y de esta forma, reformular el problema en una forma algebraica. Esto usualmente resulta en problemas de [álgebra lineal](https://es.wikipedia.org/wiki/%C3%81lgebra_lineal) de gran escala. Dos de las técnicas principales para reformular las [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) a una forma algebraica son, los [métodos de diferencias finitas](https://es.wikipedia.org/wiki/M%C3%A9todo_de_las_diferencias_finitas) (MDF), donde las [derivadas](https://es.wikipedia.org/wiki/Derivada) del problema son aproximadas por medio de la fórmula de [diferencias finitas](https://es.wikipedia.org/wiki/Diferencia_finita); y los [métodos de los elementos finitos](https://es.wikipedia.org/wiki/M%C3%A9todo_de_los_elementos_finitos) (MEF), en donde la función incógnita se escribe como combinación lineal de funciones de base simple que pueden ser derivadas e integradas fácilmente. En muchos casos, estos métodos numéricos van a ser las únicas herramientas que podamos utilizar para resolver a las [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales).\n", "\n", "### Solución de ecuaciones en derivadas parciales básicas \n", "\n", "Como resolver este tipo de ecuaciones es una tarea realmente complicada, vamos a empezar por resolver analíticamente las más fáciles para ganar confianza y así luego poder pasar a ecuaciones más complicadas. \n", "\n", "La más simple de las [Ecuaciones en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) que nos podemos encontrar es:\n", "\n", "$u_x = 0$, en donde $u = u(x, y)$.\n", "\n", "Esta ecuación nos dice que la [derivada parcial](https://es.wikipedia.org/wiki/Derivada_parcial) de $u$ con respecto a $x$ es cero, lo que significa que $u$ no depende de $x$. Por lo tanto, la solución de esta ecuación va a ser $u=f(y)$, en donde $f$ es una función arbitraria de una variable. Por ejemplo, $u = y^2 - y$ podría ser una posible solución.\n", "\n", "Subiendo un poco más la complejidad, podemos pasar a una [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) de segundo orden, como la siguiente:\n", "\n", "$u_{xx} = 0$\n", "\n", "En este caso, podemos integrar una vez $u_{xx}$ para obtener $u_x(x, y) = f(y)$. Si volvemos a integrar este resultado, podemos arribar a la solución final $u(x, y) = f(y)x + g(y)$, en donde $f$ y $g$ son dos funciones arbitrarias.\n", "\n", "Por último, si quisiéramos resolver la siguiente [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales):\n", "\n", "$u_{xy} = 0$\n", "\n", "Primero integramos con respecto a $x$ tomando a $y$ como fija, de esta forma obtenemos $u_y(x, y) = f(y)$. Luego podemos integrar con respecto a $y$ tomando a $x$ como fija, y llegamos a la solución:\n", "\n", "$u(x, y) = F(y) + g(x)$, en donde $F' = f$\n", "\n", "Como podemos ver de estos ejemplos, las soluciones analíticas de las [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) dependen de **funciones arbitrarias** (en lugar de *constantes* arbitrarias como era el caso de las [EDO](https://relopezbriega.github.io/blog/2016/01/10/ecuaciones-diferenciales-con-python/)). Por lo tanto vamos a necesitar condiciones auxiliares para poder determinar una única solución.\n", "\n", "### La condición inicial y la condición de frontera\n", "\n", "Al igual que nos pasaba cuando vimos las [ecuaciones diferenciales ordinarias](https://relopezbriega.github.io/blog/2016/01/10/ecuaciones-diferenciales-con-python/); las [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) pueden tener muchas soluciones, pero a nosotros nos va a interesar encontrar la solución para un caso particular; para lograr esto, debemos imponer unas condiciones auxiliares al problema original. Estas condiciones van a estar motivadas por la [Física](https://es.wikipedia.org/wiki/F%C3%ADsica) del problema que estemos analizando y pueden llevar a ser de dos tipos diferentes: [condiciones iniciales](https://en.wikipedia.org/wiki/Initial_condition) y [condiciones de frontera](https://es.wikipedia.org/wiki/Problema_de_condici%C3%B3n_de_frontera).\n", "\n", "La [condición inicial](https://en.wikipedia.org/wiki/Initial_condition) va a establecer el estado del problema al momento de tiempo cero, $t_0$. Por ejemplo para el problema de *[difusión](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_de_difusi%C3%B3n)*, la [condición inicial](https://en.wikipedia.org/wiki/Initial_condition) va a ser:\n", "\n", "$$u(x, t_0) = \\phi(x)$$\n", "\n", "donde $\\phi(x)= \\phi(x, y, z)$ es una función que puede representar el estado de concentración inicial. Para el problema del *[flujo del calor](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_del_calor)*, $\\phi(x)$ va a representar la temperatura inicial.\n", "\n", "La [condición de frontera](https://es.wikipedia.org/wiki/Problema_de_condici%C3%B3n_de_frontera) nos va a delimitar el *dominio* en el que nuestra [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) es válida. Así por ejemplo, volviendo al problema de *[difusión](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_de_difusi%C3%B3n)*, el *dominio* en el que nuestra [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) es válida, puede estar delimitado por la superficie del objeto que contiene al líquido. Existen varios tipos de [condiciones de frontera](https://es.wikipedia.org/wiki/Problema_de_condici%C3%B3n_de_frontera), de las cuales las más importantes son:\n", "\n", "* La **[condición de frontera de Dirichlet](https://es.wikipedia.org/wiki/Condici%C3%B3n_de_frontera_de_Dirichlet)**, en dónde los valores válidos de la función incógnita $u$ son especificados.\n", "\n", "* La **[condición de frontera de Neumann](https://es.wikipedia.org/wiki/Condici%C3%B3n_de_frontera_de_Neumann)**, en donde los valores válidos especificados son dados para alguna de las derivadas de $u$.\n", "\n", "* La **[condición de frontera de Robin](https://es.wikipedia.org/wiki/Condici%C3%B3n_de_frontera_de_Robin)**, en donde los valores válidos son especificados por una [combinación lineal](https://es.wikipedia.org/wiki/Combinaci%C3%B3n_lineal) de una función y las derivadas de $u$.\n", "\n", "### Interpretación geométrica de EDP de primer orden \n", "\n", "Las [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) de primer orden poseen una interpretación geométrica la cual nos puede facilitar alcanzar una solución general para ellas. \n", "\n", "#### Coeficientes constantes\n", "\n", "Tomemos la siguiente ecuación de coeficientes constantes:\n", "\n", "$au_x + bu_y = 0$, en donde $a$ y $b$ son *constantes* y ambas no pueden ser cero.\n", "\n", "En esta ecuación la cantidad $au_x + bu_y$ es la derivada direccional de $u$ en la dirección del [vector](https://es.wikipedia.org/wiki/Vector) $V = (a, b) = ai + bj$. Como esta cantidad tiene que ser cero, esto significa que $u (x, y)$ debe ser constante en la dirección de $V$. El [vector](https://es.wikipedia.org/wiki/Vector) $(b, -a)$ es ortogonal a $V$, por lo tanto, las líneas paralelas a $V$ (ver el gráfico más abajo) tienen las ecuaciones $bx - ay$ constantes. (Se las llama las *líneas características*.) Entonces, la solución es constante en cada una de esas líneas. Por lo tanto, $u (x, y)$ depende de $bx - ay$ solamente. De esta forma podemos llegar a la solución general para este tipo de ecuaciones, que va a ser:\n", "\n", "$u(x, y) = f(bx - ay)$, en donde $f$ es una función arbitraria de una variable.\n", "\n", "\n", "\"EDP\n", "\n", "Entonces, si por ejemplo, quisiéramos resolver la [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales), $4u_x - 3u_y= 0$, con la [condición de frontera](https://es.wikipedia.org/wiki/Problema_de_condici%C3%B3n_de_frontera) que $u(0, y) = y^3$. Podemos aplicar la solución general que obtuvimos arriba y llegar al resultado $u(x, y) = f(-3x - 4y)$. Ahora, solo nos faltaría aplicar la condición para poder determinar cual es la función arbitraria $f$. Si sustituimos $x=0$ en nuestra solución, obtenemos $y^3 = f(-4y)$. Si decimos que $z = -4y$, entonces nuestra función es $f(z) = -z^3 / 64$. Por lo tanto la solución de nuestra [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) que satisface la [condición de frontera](https://es.wikipedia.org/wiki/Problema_de_condici%C3%B3n_de_frontera) es $u(x, y) = (3x + 4y)^3 / 64$. \n", "\n", "#### Coeficientes variables\n", "\n", "Consideremos ahora la siguiente [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) de coeficientes variables:\n", "\n", "$u_x + yu_y = 0$\n", "\n", "Esta ecuación es similar a la que vimos anteriormente, con la diferencia de que ahora tenemos al coeficiente variable $y$. Utilizando la misma intuición geométrica que usamos antes, podemos ver que aquí también la derivada direccional de $u$ en el vector $v=(1, y)$ es constante, pero esta vez el vector no es constante, sino que es variable con $y$. Las curvas que tienen a $v$ como su vector tangente, tienen pendiente $y/1$, es decir:\n", "\n", "$\\frac{dy}{dx} = \\frac{y}{1}$.\n", "\n", "Podemos resolver esta ecuación como una [EDO](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) y así obtener:\n", "\n", "$y = Ce^x$, o lo que es lo mismo $e^{-x}y = C$.\n", "\n", "La solución de nuestra [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) entonces va a ser constante en estas *curvas características* (ver gráfico); y van a responder a la siguiente solución general:\n", "\n", "$u(x, y) = f(e^{-x}y)$, en donde $f$ es una función arbitraria.\n", "\n", "\"EDP" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Resolviendo ecuaciones en derivadas parciales con Python\n", "\n", "Es tiempo de nuevamente recurrir a nuestros queridos paquetes científicos de [Python](https://www.python.org/), [NumPy](https://www.numpy.org/), [Matplotlib](https://matplotlib.org/), [SymPy](https://www.sympy.org/es/) y [SciPy](https://www.scipy.org/) para ayudarnos a resolver las [Ecuaciones en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales). Así como en el caso de las [Ecuaciones diferenciales ordinarias](https://relopezbriega.github.io/blog/2016/01/10/ecuaciones-diferenciales-con-python/) vimos que existía dentro del paquete [SymPy](https://www.sympy.org/es/), el solucionador genérico `sympy.dsolve`; para el caso de las [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales), vamos a tener al solucionador [`sympy.pdsolve`](https://docs.sympy.org/dev/modules/solvers/pde.html). Aunque en este caso, es mucho más limitado que su versión para [EDO](https://relopezbriega.github.io/blog/2016/01/10/ecuaciones-diferenciales-con-python/); ya que solo vamos a poder resolver [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) de primer orden. Veamos como funciona:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true, "hide_input": false }, "outputs": [], "source": [ "# \n", "# importando modulos necesarios\n", "%matplotlib inline\n", "\n", "import matplotlib.pyplot as plt\n", "import matplotlib as mpl\n", "import numpy as np\n", "import sympy \n", "\n", "# imprimir con notación matemática.\n", "sympy.init_printing(use_latex='mathjax') " ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Defino las variables\n", "x, y, u, z = sympy.symbols('x y u z')\n", "f = sympy.Function('f')" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$4 \\frac{\\partial}{\\partial x} f{\\left (x,y \\right )} - 3 \\frac{\\partial}{\\partial y} f{\\left (x,y \\right )} = 0$$" ], "text/plain": [ " ∂ ∂ \n", "4⋅──(f(x, y)) - 3⋅──(f(x, y)) = 0\n", " ∂x ∂y " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Defino la EDP 4u_x - 3u_y = 0\n", "u = f(x, y)\n", "u_x = u.diff(x)\n", "u_y = u.diff(y)\n", "eq = sympy.Eq(4*u_x - 3*u_y)\n", "eq" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$f{\\left (x,y \\right )} = F{\\left (- 3 x - 4 y \\right )}$$" ], "text/plain": [ "f(x, y) = F(-3⋅x - 4⋅y)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Resuelvo la ecuación\n", "sympy.pdsolve(eq)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$y \\frac{\\partial}{\\partial y} f{\\left (x,y \\right )} + \\frac{\\partial}{\\partial x} f{\\left (x,y \\right )} = 0$$" ], "text/plain": [ " ∂ ∂ \n", "y⋅──(f(x, y)) + ──(f(x, y)) = 0\n", " ∂y ∂x " ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Defino la EDP u_x + yu_y = 0\n", "u = f(x, y)\n", "u_x = u.diff(x)\n", "u_y = u.diff(y)\n", "eq2 = sympy.Eq(u_x + y*u_y)\n", "eq2" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$f{\\left (x,y \\right )} = F{\\left (y e^{- x} \\right )}$$" ], "text/plain": [ " ⎛ -x⎞\n", "f(x, y) = F⎝y⋅ℯ ⎠" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Resuelvo la ecuación\n", "sympy.pdsolve(eq2)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "('1st_linear_variable_coeff',)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Calsificación de EDP.\n", "sympy.classify_pde(eq2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Como podemos comprobar, obtuvimos los mismos resultados utilizando [`sympy.pdsolve`](https://docs.sympy.org/dev/modules/solvers/pde.html) que en nuestro análisis manual con la interpretación geométrica. Otra limitación que vamos a tener al trabajar con [`sympy.pdsolve`](https://docs.sympy.org/dev/modules/solvers/pde.html) es que no podemos aplicar nuestras condiciones auxiliares al problema, por lo que para despejar la función arbitraria, deberíamos hacer un trabajo manual. Asimismo, [SymPy](https://www.sympy.org/es/) también nos ofrece la función `classify_pde` la cual nos ayuda a saber con que tipo de [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) estamos tratando. (Recordemos que `pdsolve` solo puede resolver [EDPs](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) de primer orden)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Separación de variables\n", "\n", "Otra forma en que nos podemos ayudar de [SymPy](https://www.sympy.org/es/) para resolver [EDPs](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales), es utilizando el [método de separación de variables](https://es.wikipedia.org/wiki/M%C3%A9todo_de_separaci%C3%B3n_de_variables). La característica esencial de esta técnica es transformar la [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) en un conjunto de [EDOs](https://relopezbriega.github.io/blog/2016/01/10/ecuaciones-diferenciales-con-python/) (las cuales podemos solucionar con la ayuda de [SymPy](https://www.sympy.org/es/)). De esta forma, la solución requerida de la [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) se expresa como\n", "un producto $u (x, y) = X (x) Y (y) \\ne 0$, o como una suma $u (x, y) = X (x) + Y (y)$, donde $X (x)$ e $Y (y)$ son funciones de $x$ e $y$, respectivamente. Muchos de los problemas significativos en [ecuaciones en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) pueden ser resueltos por este método. Para ilustrar como funciona esta técnica, veamos un ejemplo. Vamos a resolver la siguiente [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales).\n", "\n", "$y^2u_x^2 + x^2u_y^2 = (xyu)^2$, que cumple con la condición $u(x, 0) = 3e^{x^2/4}$.\n", "\n", "Podemos entonces asumir que $u(x, y) = f(x) g(y) \\ne 0$ es una *solución separable* de ella; por tanto lo reemplazamos en la ecuación para obtener:\n", "\n", "$$y^2(f'(x) g(y))^2 + x^2(f(x) g'(y))^2 = x^2 y^2 (f(x) g(y))^2$$ \n", "\n", "lo que es equivalente a decir;\n", "\n", "$$\\frac{1}{x^2}\\left(\\frac{f'(x)}{f(x)}\\right)^2 + \\frac{1}{y^2}\\left(\\frac{g'(y)}{g(y)}\\right)^2 = 1$$\n", "\n", "Luego, ayudándonos de la *constante de separación* $\\lambda^2$, podemos separar a esta ecuación en dos [EDOs](https://relopezbriega.github.io/blog/2016/01/10/ecuaciones-diferenciales-con-python/), del siguiente modo; primero igualamos la ecuación anterior a $\\lambda^2$:\n", "\n", "$$\\frac{1}{x^2}\\left(\\frac{f'(x)}{f(x)}\\right)^2 = 1 - \\frac{1}{y^2}\\left(\\frac{g'(y)}{g(y)}\\right)^2 = \\lambda^2$$\n", "\n", "y luego separamos ambas ecuaciones para obtener:\n", "\n", "$$\\frac{1}{x}\\frac{f'(x)}{f(x)} = \\lambda \\\\ \\frac{g'(x)}{yg(y)}= \\sqrt{1 - \\lambda^2}$$\n", "\n", "Ahora podemos utilizar `sympy.dsolve` para resolver ambas [EDOs](https://relopezbriega.github.io/blog/2016/01/10/ecuaciones-diferenciales-con-python/):" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$- z + \\frac{\\frac{d}{d x} f{\\left (x \\right )}}{x f{\\left (x \\right )}} = 0$$" ], "text/plain": [ " d \n", " ──(f(x)) \n", " dx \n", "-z + ──────── = 0\n", " x⋅f(x) " ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# EDO n° 1\n", "edo1 = sympy.Eq((1 / x) * (f(x).diff(x)/f(x)) - z)\n", "edo1" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$f{\\left (x \\right )} = C_{1} e^{\\frac{x^{2} z}{2}}$$" ], "text/plain": [ " 2 \n", " x ⋅z\n", " ────\n", " 2 \n", "f(x) = C₁⋅ℯ " ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Resolviendo EDO n° 1\n", "sympy.dsolve(edo1)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$- \\sqrt{- z^{2} + 1} + \\frac{\\frac{d}{d y} f{\\left (y \\right )}}{y f{\\left (y \\right )}} = 0$$" ], "text/plain": [ " d \n", " __________ ──(f(y)) \n", " ╱ 2 dy \n", "- ╲╱ - z + 1 + ──────── = 0\n", " y⋅f(y) " ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# EDO n° 2\n", "edo2 = sympy.Eq((f(y).diff(y)) / (y*f(y)) - sympy.sqrt(1 - z**2))\n", "edo2" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$\\left [ f{\\left (y \\right )} = C_{1} e^{- \\frac{1}{2} \\sqrt{y^{4} \\left(- z^{2} + 1\\right)}}, \\quad f{\\left (y \\right )} = C_{1} e^{\\frac{1}{2} \\sqrt{y^{4} \\left(- z^{2} + 1\\right)}}\\right ]$$" ], "text/plain": [ "⎡ _______________ _______________⎤\n", "⎢ ╱ 4 ⎛ 2 ⎞ ╱ 4 ⎛ 2 ⎞ ⎥\n", "⎢ -╲╱ y ⋅⎝- z + 1⎠ ╲╱ y ⋅⎝- z + 1⎠ ⎥\n", "⎢ ──────────────────── ──────────────────⎥\n", "⎢ 2 2 ⎥\n", "⎣f(y) = C₁⋅ℯ , f(y) = C₁⋅ℯ ⎦" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Resolviendo EDO n° 2\n", "sympy.dsolve(edo2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Entonces ahora podemos utilizar estos resultados para armar la solución final a nuestra [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) original, el cual va a ser:\n", "\n", "$$u(x, y) = C e^{\\frac{\\lambda}{2}x^2 + \\frac{1}{2}y^2\\sqrt{1 - \\lambda^2}}$$\n", "\n", "En donde $C = C_1 C_2$ es una *constante arbitraria*.\n", "\n", "Por último, utilizando la condición $u(x, 0) = 3e^{x^2/4}$, podemos despejar tanto a $C$ ($C=3$) como a $\\lambda$ ($\\lambda = 1/2$) y arribar a la solución final:\n", "\n", "$$u(x, y) = 3 e^{\\frac{1}{4}\\left(x^2 + y^2\\sqrt{3}\\right)}$$\n", "\n", "Si bien debemos realizar un trabajo manual previo, aun así [SymPy](https://www.sympy.org/es/) sigue siendo de gran ayuda para facilitarnos llegar a la solución final.\n", "\n", "### Métodos numéricos - Método de Elementos Finitos\n", "\n", "Por último, para cerrar este artículo, veamos como podemos aplicar el [métodos de los elementos finitos](https://es.wikipedia.org/wiki/M%C3%A9todo_de_los_elementos_finitos) (MEF) con [Python](https://www.python.org/). Para esto nos vamos ayudar de la librería [FEniCS](https://fenicsproject.org/), la cual es un *framework* para resolver numéricamente problemas generales de [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) utilizando el [métodos de los elementos finitos](https://es.wikipedia.org/wiki/M%C3%A9todo_de_los_elementos_finitos). Para instalar esta librería en [Ubuntu](https://www.ubuntulinux.org/), pueden utilizar los siguientes comandos:\n", "\n", "```\n", "sudo add-apt-repository ppa:fenics-packages/fenics\n", "sudo apt-get update\n", "sudo apt-get install fenics\n", "```\n", "Deben tener en cuenta que por ahora solo funciona con [Python 2](https://www.python.org/downloads/release/python-2711/). La interfaz principal que vamos a utilizar para trabajar con este *framework* nos la proporcionan las librerías `dolfin` y `mshr`; las cuales debemos importar para poder trabajar con el. Una vez importadas, podemos configurar algunos de sus parámetros para lograr el comportamiento deseado." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# importando modulos de fenics\n", "import dolfin\n", "import mshr\n", "\n", "dolfin.parameters[\"reorder_dofs_serial\"] = False\n", "dolfin.parameters[\"allow_extrapolation\"] = True" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "El problema que vamos a resolver con la ayuda de [FEniCS](https://fenicsproject.org/), va a ser la siguiente [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales):\n", "\n", "$$u_{xx} + u_{yy} = 0$$\n", "\n", "Con las siguientes [condiciones de frontera](https://es.wikipedia.org/wiki/Problema_de_condici%C3%B3n_de_frontera):\n", "\n", "$$u(x=0) = 3 ; \\ u(x=1)=-1 ; \\ u(y=0) = -5 ; \\ u(y=1) = 5$$ \n", "\n", "El primer paso en la solución de una [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) utilizando el [métodos de los elementos finitos](https://es.wikipedia.org/wiki/M%C3%A9todo_de_los_elementos_finitos), es definir una *[malla](https://es.wikipedia.org/wiki/Malla_poligonal)* que describa la discretización del *dominio* del problema. Para este caso, vamos a utilizar la función `RectangleMesh` que nos ofrece [FEniCS](https://fenicsproject.org/)." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Discretizando el problema\n", "N1 = N2 = 75\n", "mesh = dolfin.RectangleMesh(dolfin.Point(0, 0), dolfin.Point(1, 1), N1, N2)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "El siguiente paso es definir una representación del [espacio funcional](https://es.wikipedia.org/wiki/Espacio_funcional) para las funciones de ensayo y prueba. Para esto vamos a utilizar la clase `FunctionSpace`. El constructor de esta clase tiene al menos tres argumentos: un objeto de *[malla](https://es.wikipedia.org/wiki/Malla_poligonal)*, el nombre del tipo de función base, y el grado de la función base. En este caso, vamos a utilizar la función de [Lagrange](https://es.wikipedia.org/wiki/Interpolaci%C3%B3n_polin%C3%B3mica_de_Lagrange)." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "DEBUG:FFC:Reusing form from cache.\n" ] } ], "source": [ "# Funciones bases\n", "V = dolfin.FunctionSpace(mesh, 'Lagrange', 1)\n", "u = dolfin.TrialFunction(V)\n", "v = dolfin.TestFunction(V)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "Ahora debemos definir a nuestra [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) en su [formulación débil](https://es.wikipedia.org/wiki/Formulaci%C3%B3n_d%C3%A9bil_de_una_ecuaci%C3%B3n_diferencial) equivalente para poder tratarla como un problema de [álgebra lineal](https://es.wikipedia.org/wiki/%C3%81lgebra_lineal) que podamos resolver con el [MEF](https://es.wikipedia.org/wiki/M%C3%A9todo_de_los_elementos_finitos)." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Formulación debil de la EDP\n", "a = dolfin.inner(dolfin.nabla_grad(u), dolfin.nabla_grad(v)) * dolfin.dx\n", "f = dolfin.Constant(0.0)\n", "L = f * v * dolfin.dx" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Por último, solo nos falta definir las [condiciones de frontera](https://es.wikipedia.org/wiki/Problema_de_condici%C3%B3n_de_frontera)." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Defino condiciones de frontera\n", "def u0_top_boundary(x, on_boundary):\n", " return on_boundary and abs(x[1]-1) < 1e-8\n", "\n", "def u0_bottom_boundary(x, on_boundary):\n", " return on_boundary and abs(x[1]) < 1e-8\n", "\n", "def u0_left_boundary(x, on_boundary):\n", " return on_boundary and abs(x[0]) < 1e-8\n", "\n", "def u0_right_boundary(x, on_boundary):\n", " return on_boundary and abs(x[0]-1) < 1e-8" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Definiendo condiciones de frontera de Dirichlet\n", "bc_t = dolfin.DirichletBC(V, dolfin.Constant(5), u0_top_boundary)\n", "bc_b = dolfin.DirichletBC(V, dolfin.Constant(-5), u0_bottom_boundary)\n", "bc_l = dolfin.DirichletBC(V, dolfin.Constant(3), u0_left_boundary)\n", "bc_r = dolfin.DirichletBC(V, dolfin.Constant(-1), u0_right_boundary)\n", "\n", "# Lista de condiciones de frontera\n", "bcs = [bc_t, bc_b, bc_r, bc_l]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Con esta especificación de las [condiciones de frontera](https://es.wikipedia.org/wiki/Problema_de_condici%C3%B3n_de_frontera), ya estamos listos para resolver nuestra [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) utilizando la función `dolfin.solve`. El vector resultante, luego lo podemos convertir a una matriz de [NumPy](https://www.numpy.org/) y utilizarla para graficar la solución con [Matplotlib](https://matplotlib.org/)." ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "DEBUG:FFC:Reusing form from cache.\n", "DEBUG:FFC:Reusing form from cache.\n" ] } ], "source": [ "# Resolviendo la EDP\n", "u_sol = dolfin.Function(V)\n", "dolfin.solve(a == L, u_sol, bcs)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi0AAAGsCAYAAAAR7ZeSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuYbVdZ5/vfW6uq9iW3DYEQSSLQISAohiAJHBDBE2ki\nLdLSfRqDyAN4we4nPJ4+4oXQRzzt8YKnVbRzBEQu2k0bn4NogxLSUS4JF2MiIYBJMAmXTgIEQ0LI\nZSd776r3/LEqsXbN903GqLnmmnOu+n6ep56kZo055lhzV60aNeZvvtPcXQAAAEO31PcAAAAASjBp\nAQAAo8CkBQAAjAKTFgAAMApMWgAAwCgwaQEAAKPQ66TFzN5uZjeb2WceoM3vmtm1ZnalmZ02z/EB\nALCIzOyLZvZpM7vCzP627/GU6nul5R2Szsq+aGbPl/RYdz9F0k9KetO8BgYAwAJzSc9x99Pc/Yy+\nB1Oq10mLu18i6bYHaPKDkv5wo+2lkvaZ2SPmMTYAABac9T2AWn2vtDyYEyTdsOnzGyWd2NNYAABY\nFC7pr8zscjP7ib4HU2q57wEU2DoT5LkDAAC080x3/4qZPVzSRWZ2zcbVj0Eb+qTlJkknbfr8xI1t\nhzEzJjIAgLlx904vrcz699rW8br7Vzb++49m9meSzpDEpKWl90o6R9L5ZvZ0Sd9w95ujhle9/AWH\nff477/x096NbMJfrG3qq9vU9jIXAuZwdzuVscB7r/bubrgy3n3rCfM7jq/SomfTzFn3psM/NbK+k\nibvfYWZHSPrnkv6vmRysY71OWszsjyU9W9LDzOwGSa+XtCJJ7v4Wd3+/mT3fzK6TdJekV/Q3WgAA\n5mcyq7Wc5prNIyT9mZlJ03nAu9z9f8zoaJ3qddLi7mcXtDlnHmMBAGAncPcvSHpy3+PYjqFfHsIc\nPVK7+x7CwuBczg7ncjY4j+MzsRkttSxQ6nNhJy3PP/7Ixrb/effBsO3tB9cb2w56/K98YL25fS1p\nG1mr+Oap6XcWTtaeuRyn5hy0NbPl1UrzOpeLaOsb9eO0t6eRPLh5fn+1/QX2UB1RNd7oeFGNjGxc\n0bGytqtLze3RNknaE3R8zEpzZEcdG3/fPOzbHtrY9tgXPTNs+/lw6/z09f41ZAs7aQEAYMxmttKy\nQIZeXA4AAEASKy0AAAwSl4eamLQAADBAXB5qWthJy9EnHNXY9q033RG2vfXAWmPbnYea4VwpDuLG\n4dx4XFG4Nj5S+yDuPAOvNeYZMOaHvs5Q/7Lr6t+x7fXxWYyr9JzXBF5r+sj2j9quBiesJlwbhWgl\naXVl0ti26+hdYdtdR682th35iCMa2/Y9Nn627sNPO6U5ru9+YdgWw7OwkxYAAMZsqH9E9IlJCwAA\nA8RKcRN3DwEAgFFgpQUAgAFiVaFpYSctUVhr7/64SunSbfc0tq0uNcO5krR/rRmbPRAkaaNwrhSH\nUGtCu23DtfOushtZ1/yWPHfaD/08l5PHVA122sfsj9U2BFvTb01gti5cWx6knSw3f6Imq80QrSQt\n727+elneE//Kid6v9zwkfuzAnmObodsjT3h4Y9u+x50UH+txpzW2re07IWyru+LN88Lloaad9p4O\nAABGamFXWgAAGDPuHmpi0gIAwABxeaiJSQsAAAPESkvTwk5aVoIQ2OqRK2FbX28maY++/UDYNgqn\nlYZzJWnNm/vPIohbGrCNjt+lIQR/I2P6C2Zsgde437b7dxWOnd8Yairatg3XRoFZSVoKOllKgrST\noEptFKSN3mslaSV4v109In4Pjqrf7j72mLDt3uP2NbYdeVKz+u3KiY8N97fjvrWx7dDqkWFb3XVv\nvB29WdhJCwAAYzamP67mhUkLAAADxOWhJm55BgAAo8BKCwAAA8TloaaFnbREVRqjCo2StHagWf12\nPUm82v5DzWMdbLbLK+K2a5tXxC0P+MbH6iYwO+/gb6mhLrt28SbV1WvtKsQ6z2O1DczWjGEp6deC\njqP3r6idFIdroxCtFIduJyvxgntp6Darcrt6RLPK7erRcVXyXfuOamzbfezRYdu9xx/bHMMjmuHa\n5UfEFXHXdjcDvvccSu6c6NlQ36f6xOUhAAAwCgu70gIAwJix0tLEpAUAgAEi09LE5SEAADAKC7vS\nUhpuk+KArlekWKMg3EpSEvdgUH03WwKMq+eWj6suiFvetov9563vZddFqDzbdgzdVblNAquFx6sJ\nvNok/rsv6iPaX4pDszXvXxZU6U7f64LQbBbaLQ3irhzRrGYrSStHH9HYtnrU3rDtrn3NirRZRdzJ\nsd/S3PaQ4xrbfE8c5F3f1TzWoeRmiL71/T41RAs7aQEAYMy4PNTE5SEAAHYgM5uY2RVm9r6+x1KK\nlRYAAAZoDpeHflrSVZKahXIGipUWAAAGaGI2k4+ImZ0o6fmS/kBRddKBWtiVliggl1V+jEJrXhHM\nWgsDenEQ14Lqu5OkGmMUuq2pMBvvn7Vu1+8imOf141n8BTXUKrORKASbBV7D/Svatg3H2lK8fxzE\nzarcNvvoKogbva/NIogbBWwnu5tVbleOiKvcrh7dDN2uHtUM50rSroc2Q7dLxzQr30rSJNi+dPRD\nG9vWVpuBW0ny5ebrWjsQNu1dxystvy3pZyXFieWBYqUFAIAdxMx+QNLX3P0KjWiVRVrglRYAAMZs\nuyucn7n3Ln32wF0P1OQZkn7QzJ4vabeko83sj9z9Zds64BwxaQEAYIBKawttderuI3Xq7n+6PHb+\nXf942Nfd/VxJ50qSmT1b0mvGMGGRuDwEAMBON5qg4o5aackCemEQN0msrrcs/RqGAYNwriRNgmOl\n4wrDsV1V1B3VJdBiXYXeugix1ph34LWLMWTh2NL9szHUBGajyrNLyTlYWo36Td5/giBteCNB0KcU\nV/TOq383g7RRuHbathlYXT5id2PbSrBNklaDirgrR8ZB3KWj9jW3JUHcaPv6rma/vhoHhA8Ff6uv\nrcfvwX2r+bnZLnf/iKSPdH6gGdlRkxYAAMai5o+InYLLQwAAYBRYaQEAYIBqLsHuFExaAAAYoHlk\nWsZmYSct0bXAmoBeFKSTpMla2cw3+2ZbO9CsfpuNKwr9ZkFcX2v2G+0fx/PifleStjtN2zeOttel\nuwq8RvIQavsgbKPPGVSTDfcPArPTPqLAa3mV2yjcupQcqyZIG/UbvYaacO3SavzWPonCtUkQNwrY\nLu+Ngrhx4HX5yGZF2qUj4kfcLB31kOa2vcnjcKLQ7XJzXD6JX1f0Ftry/grM0cJOWgAAGDOCuE1M\nWgAAGKCaFc6dgjMCAABGgZUWAAAGiMtDTTtq0pIF7KIwnVdUz414xTebr2XhxyCIO2kGbiVpPegj\negVZkLetKAjclXnfBth7ELeDEKxU97q6CMdmFWLrxtWu3+hnP/sZL62oK8VVbvPquc3jTVaab801\n4dpo/2nbZjg1D+I2A7ZRRdzJnr3h/lHoNgrcSnHodumIo8O26yvN1+vLzW2axLcSrK033wPjiuL9\n4+6hph01aQEAYCyo09LEGQEAAKPASgsAAANEpqVpYSctWYGpsG1UdCq47ilVZFoqsiPZk6MtzInE\ni2PR8cJ+KyrG+XpNTmWYi3azuGWwfXG4bp6cHImfkFzzhOX2hexqcibF+1cUcWv7lOcojzLtNzhW\nln8Jvu8mSSZlKcyvNH9Qa3Iq2bHiTEuQB1GcX1na3cyvLO2Jn9xse5uZFAv2lyQL8itRwbh0+3Lz\nda1nmZbgfXWYiZa6n92dYpi/aQAAALZY2JUWAADGrGaFcqdg0gIAwABxy3MT0zgAADAKO2qlJQs5\nLk2iGFYcsFvTWrNlEMZbO9BsN+01eqJrHHitKVAXhW7j4nLjD9fOQheB1/RYcw7CbrfdAx2/qo/g\n9XYRmM36aPuU5ygYm+2fFXyLgrhRuFZKCskF25b3xEXgwrZJuDYM7SbF5aKAre0u25btn4Z2V5uF\n7KIicpLConHRE52jInKStDbQQnIRVlqadtSkBQCAsSDT0sQZAQAAo8BKCwAAA8TloSYmLQAADFBW\nVHEn21GTliz4FwVes+tm0ROZo2hrVjl3PQjCWvKU55qqulmYtyGp9ulJaG2IZlElsu214rZPIs50\nEcStqSab9hGGY9tXqQ33b/lE6JogbTjWZH+blD2NubbfKMwb9ZsFZidR9dys7a5muNV2xZVnS0O3\nabg2qH6bhXY9enJzEK6VJA+q32opuBkieUuLcrhDffvjgYlNnBEAADAKO2qlBQCAseCBiU1MWgAA\nGCCCuE1cHgIAAKPASouysFMcbA0DdkHbLBY7CcJ8WeA2Cu1msjBvGzVB4Mw8/1Joe6wuQrBS3R0A\nbcOxcZ/zC8xO+yirPJu91rCabBrELQvXStLSSjOwGp6D4PhSHJjNQ7vNPrLQbhzELQ/XhuNaTcK1\nQejWsraFodsocJu19eX4WGHoNqh8m7VdD9quHYzfP6PQ7VCL5BLEber1jJjZWWZ2jZlda2Y/H3z9\nYWb2ATP7lJl91sxe3sMwAQCYu6WJzeRjKzPbbWaXbvxuvcrMfq2Hl7ctvU1azGwi6TxJZ0l6oqSz\nzewJW5qdI+kKd3+ypOdI+k0zY3UIAIBtcvd7JH3vxu/W75T0vWb23T0Pq0ifKy1nSLrO3b/o7gcl\nnS/phVvafEXS0Rv/f7Skr7v7oTmOEQCAXtiSzeQj4u53b/zvqqZPCL51Xq+rjT5XLU6QdMOmz2+U\n9LQtbd4q6YNm9mVJR0n6N3MaGwAAverygYlmtiTpk5JOlvQmd7+qs4PNUJ+TlpLo07mSPuXuzzGz\nkyVdZGanuvsdWxued8Xn7v//M44/Vqcfd2yjs7X1tRbDnZqsNr+Joiq5loRYPQrXphVLy0ORswjN\nDlHrcO0MqueWhuGqQqwV4dia49VUCy4NzNaOqzRImwZxo7ZpOLYZwMyCtKVVarNxTYLAaz6usnBt\n1jYaV1TNVpK0HASMd+0Jm0ah2yyIG4Zuo/2TY2k5qHIbjFWKw7Vh5VsprH4bBWmzd0TfRur2so9f\noss/8dHq/fryt1+9RZfd/PUHbOPu65KebGbHSLrQzJ7j7h+ex/ja6HPScpOkkzZ9fpKmqy2bPUPS\nr0iSu19vZl+Q9HhJl2/t7JzTHn/Y54v6SxwAMF+nP+NZOv0Zz7r/8zf/1hvmctzt/tH2tBMerqed\n8PD7P/+9T1+btnX3283sLyU9VdKHt3XAOeoz03K5pFPM7NFmtirpxZLeu6XNNZK+T5LM7BGaTlg+\nP9dRAgDQA5sszeSj0e/0ztx9G/+/R9JzJV0x55e3Lb2ttLj7ITM7R9KFmoaA3ubuV5vZqza+/hZJ\nvyrpHWZ2paYTrJ9z91GEhQAAGKhvkfSHG7mWJUn/xd3/uucxFen19mF3v0DSBVu2vWXT/98i6QXz\nHhcAAH3LMlNtuftnJD2lk847tqNqnqQVPINt68mVsyhIG113TL/Vgoq460n+JmgaB3mlNMw7RF1V\neZxFEDZu20EQt6KabNpHRZXZeAzdVJ6N+rXgmzmtXBvuP4vQbhB4LQznZsdqG67NtlsQQrXkWDXh\n2ig0m7aN+o3CuStxQNijirZR5VspDNdqKT5fHmxfC8K1UeVbKa9WPkRd3j00Vjtq0gIAwFhQxr+J\nMwIAAEaBlRYAAAaIlZYmJi0AAAxQV0HcMWPSojgUmX2rZAHdqGWpSRLgDAvkRencAWtb0bbtXxqz\nCMfG/ZZXkw2PVRFirRtDTZC3PMTaRZC2LshbHsSdVARpoxDsUvIzFoVrq0K/STXY0oq2luwfbZ9F\nRVzbFbQNQrdRNVsprmgbhWir2wah3bVDzffbKJwrxdVz14sKtGMImLQAADBA0R8LOx2TFgAABohM\nSxNnBAAAjAIrLQAADFD2xPKdbEdNWrLqgutRlduK6rmh5FpkVv02YiO6nFkTeI3UhGDzMfRbTTbT\nRWA2G0NNYDaSvUl2EaRNg7wtq9Rm/zaT1SCw2vJYVeHalTiwGgZso4q4WRA3CswG+0tJEDfYP2sb\nhW49GVfYNqqSK8XVb6MquUqCtMG2JIc7KlweauKMAACAUdhRKy0AAIwFKy1NTFoAABggiss1cUYA\nAMAosNKiOPgXhXOl9oHRSbC/Z89Qb2kW4dZSbUOwNWYRmK0bb1nbtoHZ2n7bVp4NjzWDIG4YPK6o\ncltTpbZqXEG/VeHaIBhqK1kQtyJIWxrETY4VhmuzIG4U2k0q4oaVbqN/h6QibhiknWRVbpvbs4q4\nUaXb6B00e1tdj/YfaGiXy0NNTFoAABggJi1NTFoAABigmhXhnYIzAgAARmFhV1qi3ET41OTEXGe4\ncywi1zZPkvc7v8JsNWPoalxtn3qc9zu/Im6R7Hy3ffJytH9U7K1m/6xtWhwuylhEReCyvFGUM8n+\nbWsyLVHRuaDtUpI9qSpkV1gwLus3LBhX8eTmsIicFOZf3OL3qvUgrFITCxxofCXE3UNNCztpAQBg\nzMi0NHFGAADAKLDSAgDAALHS0sSkBQCAASLT0rSjJi1dhVCHoO2MvKtwbHisngOztccbexG3dAwt\nQ781Bd9KC7ul48rOV2G4dtpv1LabcG1NIbqoEFzUb/Y05ppxhaHbrG1pwbfkacxR6DYN7VrznK8l\n6dpou4cF48YUuUWpHTVpAQBgLLI/DHYyJi0AAAwQmZYmzggAABgFVloAABigrlZazOwkSX8k6ThN\n6+39vrv/bicHmzEmLTPQ2TdWy+R4F4HZ2ralr2GeTz2u7bemGmxpv11Vni3tczqG8qq+XVSpDUO0\n04EF25JjRZVfk37D8bYN4mZPeY6CuOlTngv7zZ7cvLKrsS0LvIaVZ7NwbPiU56BKbjKumtCuR+NK\ngrjhE50Lt0l11XP71uHdQwcl/Xt3/5SZHSnp78zsIne/uqsDzgqTFgAABqirP4jd/auSvrrx/3ea\n2dWSHilp8JMWMi0AAOxQZvZoSadJurTfkZRhpQUAgAHq+u6hjUtD75b00+5+Z6cHmxEmLQAADNB2\nMy0XX/UFXXLVFx64b7MVSX8q6b+6+59v60A9YNKSmGeF2UxNkLY4lDnQyrNRiLZm/6ztLCrPRmqq\n0batPJuPoSxIO4uAcRdVatPwdRjETb6/onBr2yBuNq7oNWTh2ppKvUGYOOw3CMFKSeC1osptVLl2\nuj0Ix04q9g+q3EaVb6W4ym1WETfaHBW/rSmIO6ZwbonveeJj9D1PfMz9n//aez502NfNzCS9TdJV\n7v7G+Y6uHSYtAAAMUHqHXXvPlPRSSZ82sys2tr3W3T/Q1QFnhUkLAABD1NGkxd0/qpHeiDPKQQMA\ngJ2HlRYAAIaou+Jyo8WkRd2EU7sI0Va3bRmkrQqmdhR47aJCrNRNkLam8mymiyBt1f6zqFIbvd6a\nyrNdBXHDcQXh2lkcq6aibdhv8605rXIbvK68Im7Qb03AN9o/CdfG1XeTirhBEjYLx3qQsK0J0tYE\ndPtWUyl8p2AaBwAARoGVFgAAhqi7u4dGi0kLAABDxKSlgctDAABgFHbUSktXIdYw1DnncG1pYCs7\nVhSkncVrqAnSlvabhZy7qvAaKa1GW3usToK0abg2OtYMqtSWVp6dRRC3MFw7PV7wb9Y2XFsT2k0C\nr2Hl2CgEm/zb1FS5jUK3ecC3MAychWuDtlVVbuNRaT06VhjOHVHiNjGLyuqLZkdNWgAAGA0uDzUw\naQEAYIiYtDSw9gQAAEaBlRYAAAaITEsTkxb1H7ptG67N+q0J186zSm3barTpuFpWqW0bpJ1JELeL\nIG1FlduZhGOjPqJjrcQVYsNKu+mxCsO1WR/ROUzHVXOslpVjgxBsGq5tGY6NArfp8cIqt/H+67Jm\n2yReGwV084q40bGCdvHu48LloQamcQAAYBRYaQEAYIhYaWlg0gIAwADxwMQmLg8BAIBRYKUlkVZd\nbRvErahCWhqu7W5c7cK1Nf3WVH2NQrSz6Ld1ELci8BqGTdVRkHbeQdzwNQTnK6iSm7ateQ1R5dqk\nbVW4trRybdK2qqJtzbEqwrFh6LYm4BuOK35d6zXh2mhbUtE2rJ5bkbpdH1NEl7uHGpi0AAAwRGRa\nGpjGAQCAUWClBQCAAUovi+5gTFoAABgiMi0NCztpqalC2snxK4KWWQC0bUXbvsO1NeMaarh2ur0w\nhDrPcG06hnYVYtNbLFtWqa0K1wbjzcdVERCOzkMUmE0q18bB1GbVVympaJtWxA2q59ZUuQ1Dv+VV\nbtsGfLMQ7FqwPap8K8Xh2qjK7fR4ZUHaLPQb91neFv1a2EkLAABjxuWhJiYtAAAMEZOWBiYtAAAM\nEZmWBiYtqntyc3l2ZL4F48KcyGrzn3eeOZXp8coyPF3lVJZaPvF32m9ppqXdU4+z7a0zLbPIqUQ5\nk5bF9GpyKlVPbs6yI6UF25KcSvQaaorLVRWiK33yc7p/N0+Ejp7cvJYEQqLsSRYdCdsmjaOsyzqh\nlB2DSQsAAAPEs4eael17MrOzzOwaM7vWzH4+afMcM7vCzD5rZh+e8xABAOjH0mQ2H1uY2dvN7GYz\n+0wPr6qV3iYtZjaRdJ6ksyQ9UdLZZvaELW32Sfp/Jb3A3b9D0r+e+0ABAFgs79D0d+/o9Hl56AxJ\n17n7FyXJzM6X9EJJV29q8xJJf+ruN0qSu98y70ECANCLju4ecvdLzOzRnXTesT4nLSdIumHT5zdK\netqWNqdIWjGzD0k6StLvuPt/mdP4QqWh21kUjFtaaf7ztC0Olxe9C/YPjj/tt/z1ti16F56DkT1N\nOQyR1hyrZZA2P189P025JmBcVfAtC+IGAduW4dqZFKILw7Etw7VpIbvyJ0JHodso8JoVcQsLxlWF\na+O2kfgp0UnbEWV251kQdSz6nLSUfOusSHqKpDMl7ZX0CTP7G3e/ttORAQCAwelz0nKTpJM2fX6S\npqstm90g6RZ33y9pv5ldLOlUSY1Jy3lXfO7+/z/j+GP1tBOPm/mAAQA7z2Ufv0SXf+Kj8z/wNi8P\nffjST+ojl35yxoMZBit9jsPMD2y2LOlzmq6ifFnS30o6292v3tTm2zQN6z5P0i5Jl0p6sbtftaUv\nv+rlLzi8/6paIm0v2XB5iMtD+f5cHsrGxeWhtO2ILg9FzxiS4ucMpW2jfpOHD0VtD0XHqti/pq0k\nnXrCPrl7UtBnNszM1667dCZ9TR77tMZ4NzIt73P3J83kIHPS2wUzdz8k6RxJF0q6StKfuPvVZvYq\nM3vVRptrJH1A0qc1nbC8deuEBQAAlDOzP5b0cUmPM7MbzOwVfY+pVG8rLbPUdqVlElSOras82/yL\ncxZVbmtWWqJVlaqVmoonL0eVdmtWcHqvXKsZrH60rXJbsXpiy6vlbSsq9batUltVAXhS+IRlJU8t\nrlg9qaoc29HqSVgRt2VF2/C8SK1XpqIVFal8VSV7cnPYNq2eO/t+F2Kl5frLZtLX5OTTOx/vvFAR\nFwCAAUonqiNgZkuSHiPpoZreeHOzpJvd/UCbfpm0AACA1jYKwr5C0oskPVXTOcZtktY0nbxMzOzT\nkv5c0tvc/abaY4x3GgcAwCKzpdl8dD3MqZ+VdImk4yT9qqQTJa26+3Hu/i3uvkvTicvPaDr3+Esz\ne6OZ7a05FistAAAMUZaZGhAz2yPp7ZL+StJT3P1g1tbdvynpQ5I+ZGb/UdMVmfPN7Kfc/cslx9tR\nk5YorFqr9KmbVbcmJ7cWl97GnPUR7T9JblPtKiAc3t4cBUuHGq7N2q6Uv4beb0Pu6Bzk4dgghBqG\nTcvPQXptvyq0W3h7c024NjtWFLpte9t2zS3TybFKb2OWykO32a0c0U0e2X0fYfXcpN/SqrzrRfVL\nB24cFXFfI+lcd/9CzU7uvi7p3RsPQj5X0v9Rst+OmrQAAIDZcfdfbrn/LSqcsEhMWgAAGKQx3z3U\nFSYtAAAMEZOWBs4IAAAYBVZaElXP6Kl4lk7Vc3dahnaj0G1nzzmqqGjbOixac6y24Vqp9TN+6gLC\ns69Sm56DOVapbR2uzZ6703ZcHYVrq/otDDln+3v0PKGKZ/wkhWeLQ7dZ5dooSJsdqya02zZem41h\nkBZopcXMXiLpmZKu1rRGy34ze6yk75P0NXd/T0k/TFoAABiiBZm0mNnrJb1S0wcjP0XSq83see5+\nnZntl3SDCq/8MGkBAABdeoKkx7v7PZJkZk+W9Dtmdo6ktK5LhEkLAAADtEB3D11634RFktz9U2b2\nw5L+g6S/qOloYc4IAAALZSRl/At8ycxeaWY3mNl3SJK773f310k6VXktwQZWWpRVgy2rfJvuXxPk\nnUFotzR0W1V9N2sbBQeTcGx4HoOKuGlYtKZtB+FaKQ7ShudgFuHalq+3pkJs71VqOwr9hhVipfjf\nvG1F3Zbh2rSPYFtUzVaS1iuq1EYh1CxIWxq6zX7brAVfiAK3WR9Zpd5IRdPk+GNK546Pu7/HzE6W\ndI6kz2352pvN7LOlfTFpAQBgiEbw7KFS7n69pOuTr320tB8mLQAADNEwLu1si5mdKOnb3f3CWfY7\n3jMCAMACc1uayUdPfkPSBWb2rPs2mNlrzOz5bTpl0gIAAGbtM5KeJ+ny+za4+3+SNDGzH9lupwt7\neSgNL7ZUU/02UhPaXVoNqtTWBHyDIG3UZ9avBaHQaSdB2zTwWljRdgbVZGsq7RaHa7M+2p6DmnBt\n22qw2XXxvqvUZuOqCry2C9KGx8qCvIUVdaW6MLFHYeIoMJuFWIPNaeA17DdsGvYRBWbTyrWF+2d9\n1ISJoyBt23DuIHT0e2xOvibpCHffv3mju7/PzH5xu50u7KQFAIBRG3GmRdKXJL3bzNYlfUTShzf+\ne6OkR2+301GfEQAAMEgvk/RcST8p6SZJPyXpk5K+KOmD2+2UlRYAAIZo3Cst17j7ZZIuk/RuSTKz\nR0p6saQ7ttvpqM8IAAALa9wVcW8xs9O3bPuqpPdrWgV3W3bUSksWmK2pPFsaus3aRdVva8Y1WY0r\nz0ah26pKvVEl1ijYKpWHa6XyCq9J6DeuqJtU360ZV2G4VkoCtjXh2kkwhllUgw3OTV012Z6r1Gbn\nq6sqtaUuNiuFAAAgAElEQVRh4DmGa6W4om0Ujq0J19aEWGuq1Latcpuegyi0mwZ84+2lFiKgOwLu\n/vtm9hNmdqa7//rG5jMlXSjpXdvtd0dNWgAAGIuxPzDR3d+6ZdMHNS3l/7Ht9rmtSYuZvUTSMyVd\nLelt7r7fzB4r6fskfc3d37PdAQEAAI0909Lg7muSfq9NH9WTFjN7vaRXSvpbSU+R9Goze567X2dm\n+yXdILIyAABgCzP7F5JepOmtzy7p3e7e6QMTnyDp8e5+z8YAnizpd8zsHEkHt9EfAADYqsMHJprZ\nWZLeKGki6Q/c/Q2dHexw75S06u7HmNluSS81sxe6+6+U7LydScul901YJMndP2VmPyzpP0j6i230\nN3pVQd6KtktB0LImIByGc7PAaxS6zdrWBGkLK9pmod9wvBVh4rbhWil5vW3DtWn13dlXqW0dTE3a\ntq5S21W4Njm3pUHaKEQrxQHOKEA6bVsTLC1r21W4Ng/HlvVbU+U2PV9JH+G4CqvfZud7VDq6PGRm\nE0nnaRrpuEnSZWb2Xne/upMDHu4l9/3PxlziD8xsX+nO2zkjXzKzV5rZDWb2HRsH3u/ur9P0Nqbs\nexgAAPTvDEnXufsX3f2gpPMlvXAeB3b3i9z9oi3bvlG6f/VKi7u/x8xO1jQB/LktX3uzmRVfmwIA\nALEO7x46QdP86X1ulPS0rg42S0WTFjM7UdK3u/uFkuTu10u6Pmrr7h+d3fAAANihupu0jPbiWelK\ny29I+mEze7a7XyJJZvYaSVe5+/s7G12Poqc5S0nBtihPUpE9mSR5juiJzFFOJdteVTAualtRmC1r\nG2ZHKp6wHPabHSssRFeRy4lyKtPGjU1hsbRZPE057LdlsbOajMc8C75VFKebZ8G3qNibFL/Lt82p\n1PQ7z5xK1m+kpmBcXS6ncAALLMtXPZiLL75YF1988QM1uUnSSZs+P0nT1ZbBs+wb9rBGZq+VdLmk\nj25+zLSZvUDS0e6+7ep2s2BmftXLX3DYtuiX+CSYBGRtl3fvitsWTiTmeaxsO5MWMWmRmLTc37Zw\n0lIRFmXSMv9JS9sgbtg2GVk2hlNP2Cd37+7WHk1/r929f/+DNyywd8+ew8ZrZsuaxjvOlPRlTUuY\nnD2nIG6DmT3O3f+hpG3p2tPXJB2xecIiSe7+PkknV44PAAA8CPfZfDT79UOa5lIvlHSVpD/pa8Ky\n4WWlDUsvD31J0rvNbF3SRyR9eOO/N0p6dOXgAADAg8hW/WbB3S+QdEFnB9jEzN4u6TEP0ORJmpZN\neVClk5aXSXqupEdJeo6kn5L025Lu3vh/AACAyK9oOil5p6TostprSjsqnbRc4+6XSbpM0rslycwe\nKenFku4oPdhQhVmIrG1hcbjsib+lT2OuOda0beFTiyuKwFVlWqLsihTnV8Jtyf7ROUiOVVrIbrq9\noghb4bkNcy5S3dOUa7IypdmRLKfSRX5GKs6k1BSMm2fBt5qnKWda51/CdvH+Na8h3D/ZXnpua3Iq\nmbY5k0UN8i7AS5A0vePYzP7O3T8Sfd3MTintqzTTcouZnb5l21clvV/TgnIAAGCG1n02H0Pg7uc9\nwNe2Pg06VTRpcfffl/RkM/uFTZvP1PQpz8UzJAAAsPjM7EQze96s+y2uiBvMhD6oafr4YzMdEQAA\nqLrMN0Cd1Hfbdrk9d19z999z9yu32wcAAIiN/PLQZyQ9T9Mab5Ikd/9PkiZm9iPb7XQ7T3kerSzE\n2rZtVP22syc3Z+HYKNwaBFazgHAYNq14InRVkLaronelReCk4nBt2sc8w7WTlgHfWTx9uoMgbVVR\nslmEY3t+mnI+rr6PVV5MLzz+DIq4le7/QMcrVTMGtJLWdzOzX9xupztq0gIAwFiMfHrVSX03Ji0A\nAAzQUO782aZO6rsxaQEAALPWSX23zp57DQAAts/dZ/LRk07qu7HSovYB3arKtS2r52ZByTB0Gz4h\nuaIibhKODYO0bQPC2f7Rk5drAq9JteM4hFoeeK0J14Zta56mnIyruEpt24q6kjz4vmsbpG0bos36\naBukbRtsrTlWpq4ibkW4Ndy/ePeqKrOlT2NO9+9oXF3s35Wa78WhcfffN7OfMLMz3f3XNzafqelD\nGt+13X6ZtAAAMEBDnUyV6qK+G5MWAADQOXdfk/R7bfog0wIAwACNobicmZ1rZo9osf/DzOyNpe1Z\naQEAYIBGUsb/P0t6i5l9QNK7NlZTHpSZmaR/JemlqrgFemEnLVGV2hpVQdqKcG1Y5XY1/mewIPyY\nVo4trGibBV5t1+6gbXloNwrcShWh2yhwqyzEmoRro8qxWdXWoN+0em5UzbUmXFtaUVcdValNj1VR\npTb4c22tg2q0Ul1YtCZIW/r+P4tj1Ry/NEhbEzCuHUPYb2Fps6GGa6l8Ox/ufoeZvUzTrMonzey/\nS/qopEvd/fbNbc3sCElP1bRuyw9J+itJZ2+tmvtAFnbSAgDAmI3l7iF3PyTpjWb2TkmvkPSzkp5l\nZi7pdk3n3A+RNJH0N5L+UtIPuPuNtcdi0gIAwACN4+rQP3H3b2ha9fa3zWxV0vGSjtM0P/uPkr5a\ns6oSYdICAABmyt0PSPqfGx8zw6QFAIABqikcOAZm9vOSbpP09o1LStV21KQlrTwbVbTN2hYGcZeS\nSqxh9dwswNm2om1hNdpp25ZVbnftKR5XVQg2OI95YLawcq0UhlvztoVVamsq9dZUqU0r4pYFabM3\nvqpwbAdB2rTKbXj8+FhjqjybjiF8DRX7Vxytdbh1AJVn2wZsxzQPGNFQGzbuJrpR06c7f9jdb3T3\nN5jZYyT9mqa5l2o7atICAADm4h2S/o2k35T0cDP7vKYTmKskPXG7nTJpAQBggLouDNcld/8TSX8i\nSWb2REnPlvQ9ks6V9Krt9sukBQCAARrTpawH4u5XabrC8iYze5ymE5htoYw/AAADtC6fyUcfzOyh\nZvYvzeyRm7e7+z9oehv0trDSUqk0SJtVuY22p5VngxBqWmm3NHSbVa6tqXIbhG7T0G5p6DYJLseV\nZ5NgatvquW0Dvkm4Nm5bXqV2LVkjLg3HpoHXqG3Yspsg7Tyr3E6P12y8CIHZvsOxQ6g8uygrEgvm\nXZIeLekUM/trSe+W9LcbXzttu52y0gIAwAC5z+ajJx939ydIeoqkqyW9XtIVki6V9Bfb7bTXSYuZ\nnWVm15jZtRv3b2ftTjezQ2b2onmODwCAvvTxlGcz+9/M7O/NbM3MntJi+FeZ2c9IusHd/3d3P1HS\nwyTtc/d3brfT3iYtZjaRdJ6kszS9/elsM3tC0u4Nkj4gKSluAQAAZuAzmj7M8OI2nbj7n0o6X9IP\nbtp2q7vf06bfPjMtZ0i6zt2/KElmdr6kF2q6jLTZqzW9Fnb6XEcHAECP+ri04+7XSJJlBTDr+rpJ\n0h+27miTPictJ0i6YdPnN0p62uYGZnaCphOZ/1XTSUsn/4RRRdtoWyaqfpsGZsPqqllF3ChIG4d2\n47ZlVXLTti2r3EqST5pta0KwVeHaMOCbVbkNAsJJOLa4em6y/3qwQFhTpTYL4pYGabsKvLbtt6by\nbNvA7AMdr7F/crSuwrFDDcISbu3fEELOQ9PnpKXkX+ONkn7B3d2m0z4uDwEA0IKZXaTpE5i3Otfd\n3zfv8dToc9Jyk6STNn1+kqarLZt9l6TzN5apHibp+83soLu/d2tn513xufv//4zjj9UzHvPIrU0A\nAKh22ccv0eWf+Ojcj7vd1a4HG6+7P3ebQ+qd1TxcbKYHNluW9DlJZ0r6sqb3b5/t7lszLfe1f4ek\n97n7e4Kv+VUvf8Fh25Z3Ny+XTIJtkrRyRPMySNY26rdm/2i7re4O24b1ULK2wXbbFWzLaq+E+3N5\nKDteZ5eHgm1cHuLy0AMei8tDc3fqCfvk7p2u/JuZX3HjbTPp67QTH1I9XjP7kKTXuPvfzWQQM9Lb\n3UMbj6U+R9KFmpb3/RN3v9rMXmVm234uAQAA2B4z+yEzu0HS0yX9pZld0PeYNuu1Iq67XyDpgi3b\n3pK0fUXb40UVZvO2SZC2MLSb7h/81R5WrlVS5TZbKVkpa5uu1ARjyI7lk2C8NZVng/3DlYuk37BP\nqbz6btZvRUXcmtWTtWjlIG3b3NbVika0glO30lJR4TXYVrMiMotqsqWrH0Oocluj5t8B0tIM7oqZ\nl7XsB7JD7v5nkv5s/kcuQxl/AAAGiAlpE2X8AQDAKLDSAgDAAEWXlnc6Ji2K8ydLHWVawicnJ3mQ\n4ic3K35SdLgt2z/IuoTZFSnOg2RF76KcSXRHUZYnmeOdRlFORYqXaKM3k/RpzGGfYdMk/xK3LX3y\ncpSTme5ffpdOF084rntd7Y6VafuU6EzbLMKiXhboKk9SUQs0PLdDzbks6vdBG0xaAAAYoD6CuENH\npgUAAIwCKy0AAAwQl4eamLQAADBABHGbFnbSUlfwrfwqWdQ22raUFYGLngidPXk5CpFm/Yah3Yqn\nPNeEWIPQbevQbkfhWk/aRqHZ7A0ialtVbj9sGzatKq1fOq5ZFHGrK+Mfb293rG5CsFX9dvQLZOx/\nTdeEWLNzOGkbhF0v3z/6NZD9Gww1oLuTLeykBQCAMat5htVOwaQFAIABylZvdzLuHgIAAKPASgsA\nAAM09rxTF5i0JLJw7lIUpI2SXUkAVGHl2uzJzeVto4Bu+ETn5V3h7mGV2qTKbdVTnqM+2j7lOWkb\nVbTNllfjIG48hNK2s6lyWz6u0iDtLCrP1vUbVQuO2xbvX/HmTeXa7kTB1LUkUF1TpbZGFNodU5Xb\nGtnP/k7G5SEAADAKrLQAADBAO20lrwSTFgAABoi7h5q4PAQAAEZhR620LCXh2prquUurzVO2tNLc\nFlWjlZIgbRJ4jYK0lrXdtae5MQrXBtukyiq3Yb/Jt1IUuo32TyviNoO4h9JwbTNV2TZcK8WB1Zpw\nbdR2NqHd5rb2VW6zc1vRNny97cKxXVWu7eoP2Xn+hTxZ6iZwGnUbndusmu36WtRn3DYK89YEeUvD\nudMvNNt2FRpui8tDTTtq0gIAwFhw91ATkxYAAAaIlZamgS6KAQAAHI6VFgAABmidu4camLSoLogb\nVcq1qBpsUiE2CuLWhHZtV1DlNmkbBWmzwGwYup1FaDeqaBu0jarZStLB4KJuWuW2MDA77aOmbRCO\nLexTioOwWduaIG0YJu4oMNtFOLbm/bgm2Fo3rvIxDNXBinMzqcjsllaUXU/W60uDvFIS5g0Cs9OO\ny15vFhAek0X4/pw1Lg8BAIBRYKUFAIABIojbxKQFAIABqqk5tFNweQgAAIzCwq601IRrw/0rqucq\nCsxGlW+lsPptVuW2pt8wdBu1ralym4yrJrS7HmyPQpWHgmq2UhxEy4KHbcO1Wegt6qKrKrdR8Lht\nkDYNCEdtk2OFVYGTF1y6pJ1XIG7312VXod2hKg3MStJ6RfXc0qbZOYzGlVfvrfh3iAK6heHczCS5\nEaBvfdw9ZGb/j6QfkHRA0vWSXuHut899IAlWWgAAGKA1n81Hpf8h6dvd/VRJ/yDptbN+XW0waQEA\nAJIkd7/I3e9bn71U0ol9jmerhb08BADAmA3g8uUrJf1x34PYjEkLAAAD1NXdQ2Z2kaTjgy+d6+7v\n22jzOkkH3P2/dTKIbWLSImkpCNcuTZKKtmHAN2ibhFij6rdZuNZ27Qn63RW2jUK7UejWZ1HlNugj\nCtxK0qEgSBZWuU1+NqP9518Rt6xt9ldR/HrLQ6zR/vkYgnYVgdlZhGOj43VVpbbtX6JDLZO+VBGY\njUKsWTg3CrFnVXJLg7TZKYxeQvbvtRLeJNEunJsUJQ9lYfWaPrpQEyrf7NpP/o2uu+LS9Ovu/twH\n2t/MXi7p+ZLO3NYAOsSkBQCABXLKU56uU57y9Ps/v/Adv1u8r5mdJelnJT3b3e+Z/ejaYdICAMAA\nbXelpaX/LGlV0kU2XW37hLv/uz4GEmHSAgDAAPUxaXH3U+Z+0Arc8gwAAEZhR620ZBVxo+q3S6vx\nqVmKAq9hldskXBu13bU7aRtUqc2CtIWh21lUuV1bam7PwqKlQdqDSRAu3H+O4VqpvErtwayqb7A5\na9tFkDatqFsRmK3pNxIFXucZuJ0eb5ih20kUmq34CzsO3cb7R6HbrEpuaZA2rcgbvt0mbYMfkjic\nK5UGdLNTGJ7vgerp8tCg7ahJCwAAY8GkpYnLQwAAYBRYaQEAYIBYaWli0gIAwAAxaWla2ElLFK5N\n20ZVbrP9g3BsFJi1laQibmGQV0rCtS3b+nIc+vWg0u66xeUgo6qaByqqth4K2kZ9TvdvbovCuVL7\ncG0eji2rUpu+horAa9RHV0HcmnBsXUXb4PVWvPlu46m0hx9roIHbTBpk3SKrXBv2mYRr1ytCv2FF\n3GD35SQwWxXajfpIytRGAd1oXGmVW6sYFwZnYSctAACMGSstTUxaAAAYICYtTdw9BAAARmFHrbRk\nOZX4yc1Z22bOI8y0rCYF44JCcrYSP7nZw+JyWaal7CnPWcG4KL9yILkoHBWCy3ImUUYjfPJzRU4l\nPVZFTqUu0xLt32zbVaYly4OU5lfSJ0q3zJnkhejK+p3FX5GL8Jdo9OTkNu0kaRI89Vgqz6lIcS4m\nzMQo+bkJ9l9J/kyOfp7CnIukpeDfPDo1i5BTWYTv71nbUZMWAADGIvsDbSdj0gIAwACx0tJEpgUA\nAIwCKy0AAAwQKy1NO2rSkoVrJyvNcOrSSnJqouJuNU95DgK6ebi2PEirlaDfoGBc9IRmKQ7dZgXj\n5hvEbW6/51D7cG1Ncbi4uFx5ELdtEbcwpJj0URPEDfus+HeoeUPtqm2Ntv3WBGG76DdrVxXQDdpm\n+6948/1yPSjMtubl+2fiIcTf9xMLxhWMIS2qOKJ5wFCfTN4nLg8BAIBR2FErLQAAjAWXh5qYtAAA\nMEBMWpq4PAQAAEZhYVdawsq1aZXboCJu8jTlsPptVOU22CZJCsKxnlbEjdqWP6U5Ct1mVW6j0G0e\nxA3aJn8RHFqL9g9CrMmxotDtLMK1NUHaQ0HbaLhpRd2WQdzs36yLcGxNELem37Z9dtlHF9qGY2va\n1RxrdTl4QnKy/9pS89xGbfPAbfR9G7ddCR4on/3TRj+n8ZOb4/1XVB7aVdB2nob6/d2nhZ20AAAw\nZmvJH0I7GZeHAADAKLDSAgDAAHF5qIlJCwAAA8SkpWlHTVqWaoK4QZVcKa5ouxRsi9pJkgcVbdOK\nuIVVbqXy0O09SQI0Ct1mVW6jtlmI9cCh5vZ715rp3Loqt/Gx7g1SvzVVamtCuzVB3LaB1wPJuWkb\neO2qym3bJ9Mu6ht1TWB2uaJybc2xonNbE/CNgrz5b5Ho/TbLaNRUzw1+doNt+TloblshKDEaO2rS\nAgDAWLT9A2ARMWkBAGCAFnXVsQ0mLQAADBCTlqber+SZ2Vlmdo2ZXWtmPx98/UfM7Eoz+7SZfczM\nvrOPcQIAsOjM7Jc3fud+ysz+2sxO6ntMm/W60mJmE0nnSfo+STdJuszM3uvuV29q9nlJ3+Put5vZ\nWZJ+X9LTH7TvKFy7FM/RllaC05BVxF1pbg+r3yaB2ZpwrS83264nbaOwZhS6zarc1lTEDQO+SVj0\n3mhcYZXb9uHaOLRbE/AtD7xG57urcG1NvzXXwFsHcdMqomXHmkXbvtWEY2v6mFi70G1XQdyq78/l\n5vb1rHpu+Jsoeb8OQ7fl41oP2g5VTz8Lv+Hu/6ckmdmrJb1e0o/3MZBI35eHzpB0nbt/UZLM7HxJ\nL5R0/6TF3T+xqf2lkk6c5wABAOhDH5MWd79j06dHSrpl7oN4AH1PWk6QdMOmz2+U9LQHaP9jkt7f\n6YgAANjBzOxXJP2opLtVcGVjnvqetBRPI83seyW9UtIzo6+fd8Xn7v//M44/Vs8+7iGtBwcAwGUf\nv0SXf+Kjcz/udldabrv2Ct127RXp183sIknHB186193f5+6vk/Q6M/sFSb8t6RXbGkgH+p603CRp\nc8jnJE1XWw6zEb59q6Sz3P22qKNzTnt8JwMEAOxspz/jWTr9Gc+6//M3/9Yb5nJc3+akZd/JT9a+\nk598/+dfuOAdh/fr/tzCrv6bBnZ1o+9Jy+WSTjGzR0v6sqQXSzp7cwMz+1ZJ75H0Une/rrTjpSB0\nG4VzpTiIa8vlFXGjbb6ShWuD7UE4N2sbBVulOHR7b1CNtiaIm4Vr7z7YDMfeG5WZTPqoCddGrzeq\nUJv1EYWGpbogbWnbmv2zwGxV0DEIwrYNvNbtX/4E2rbX5ocQzm0buq2raNt8r6oJzGZto+/FsMqt\n4qq8bYPiiu9viCW/nSbNtw+tBWOtCdzWhMoXnZmd4u7Xbnz6Qkn5kk0Pep20uPshMztH0oWSJpLe\n5u5Xm9mrNr7+Fkm/KOkhkt5k00T9QXc/o68xAwAwD9HjRubg18zs8ZLWJF0v6d/2MYhM3ystcvcL\nJF2wZdtbNv3/j2tAt1sBADAP3sMKkLv/67kftELvxeUAAABK9L7SAgAAmrYbxF1kO2rSEla+lWRB\n9dsoXCtJtntvY5uv7GluC6rZSpKvNvdfD/aXpHuDnGMUuJXi0O09QQg1aicl4dqKIG5aETccQ3m4\nNnoNWTXZaHtN5dmattFrmEVgti7oWD6G0mO1rag7i7Zt1RxrFhVt2/YbtV1eKg85R/uvLk+K22Zh\n9dXgxoXo+2Mm/7ZBQHcpqQo8CSriLgdjWEmGNabQbU+ZlkHj8hAAABiFHbXSAgDAWHj5gtuOwaQF\nAIAB6uPuoaFj0gIAwACRaWla2ElLVP12sppUud0VVLkNtk23B6HboPptFq6NQrsHk2jR/rCabPxN\nHLWNwrFRiHa6vbxtXOU2CeIG1W9rwrX7DzT3rwniZiHBaLxZILE0HJtXxG2+hqxt2yBs+yq33bRd\n1LsgrGXgtqZtTUXc7OcxqnKbVcRdC8K8bUPle1bjgHAkfb3B29JK0PZgsv968HrXWdEYjYWdtAAA\nMGaLOtlvg0kLAAADxKSliVueAQDAKLDSAgDAAJG1aVrYSUsUxM0r4jYDuku7j4g73tXc7itBldyg\n8q0kHbTmGPYHIdhs+71JRdy7DjTb3nngUGNbHsRtbo/CuVJ5uFYqD9JG7bK2WcgwCtLm4diobXnA\nt6YyaF2V2/K2pUvH2R0IbW+nrFm6XtS7IJYqwrWHKtpaUA02O1YUBs5CrFHoNvt52rXc3B7tn1Xf\nXVsvD91G0iDuriB0G3x/7Uq+v+dZnbktLg81cXkIAACMwsKutAAAMGastDQxaQEAYIAW9bJqGws7\nabGlINOymmRagvyK7YkzLb4aFJcLtq0nmZb9QfbkrizTEhSSuyvJftwR5FfuDNrWZVritlH+JMuk\nRNvDwm4VOZW8uFz5k6rrisOVZWWyv4qiN54sT9I2J1LzrJKaTEvbN89F/YuxprhcTf4lyrRYcjE/\n6jfLzxw61Owkew0HwvxKc9ve1Yp8VsX3XE0xvbC43Hp8wqL8SxLJwwAt7KQFAIAx49lDTUxaAAAY\nIJ7y3MSkBQCAASLT0sQtzwAAYBQWdqUlCt0urSZPbg5Ct7b7yLDt+mqz7fquoxrbssJs8dOU49l0\nFK69/Z7mNikuJFcTxL3z3ub+NeHau7O2wfGiYGsUopXqQrtR8C9rGwVDs79q1oOUXk0Itm3gtSbE\nGr6uiuPPIjC7qKHbSE0Qtyq0GwZxy4vLZaHftUnzmzRtG4R2o3BuFmBfDYrTzUIUxI0CwruSQpzR\nz8NQK8/upJ+lUgs7aQEAYMyYtDRxeQgAAIwCKy0AAAzQUC9b9YlJCwAAA8TloaaFnbRMgic62644\niLu0txmkXQ+e5jzd3gzo3rPeDIbdGVS+zbbfuv9g2Pb2IBz7zXvitncEQdi7onBuFuQNtteEa/cH\nx5LKg7RZYPZQsL1tYHbatjyIG71xdBV4bV0Rt/WTm1vtLmln3aZZVeW24mJ826c8Z6HdKOC7FIRY\nJWkyCb7HoyepJz+7q6vNpzznTzcv/1VUGsTdPUkq4q4FYWJWNEZjYSctAACMGSstTQRxAQAYoPV1\nn8nHdpjZz5jZupk9dMYvqxUmLQAA4H5mdpKk50r6Ut9j2YpJCwAAA+TuM/nYht+S9HMzfjkzsbCZ\nlrAibhC4laSlI/c1th3ac0zY9sDynsa2b97TDKZ+Mwmx3np3M0ibBXFvC0K3dwThXEm6Peg3Ctdm\n+0dVbrNwbRTQzYK0B4K20Q9RFuarCcyuR6Hd5Ae2NFybHa/mjSAKt9Ys2c4itBtpezvlEK63R2Oo\nqTzbla6q35YeKwvtLk2CtgeTtkG4dRKEWydJkDcK7UbBeikP6EYmwfmKgrh7giCwJB0MQ/TFh5+r\nPn7GzOyFkm50909HgfC+LeykBQAANJnZRZKOD770OkmvlfTPNzefy6AKMWkBAGCAthui3f/lv9c9\nX/n79Ovu/txou5l9h6THSLpyY5XlREl/Z2ZnuPvXtjWYGWPSAgDAAPl6HDN4MLuP/zbtPv7b7v/8\n9k/+f2XHc/+spEfc97mZfUHSd7n7rdsaSAeYtAAAMEDbnbTMcgh9D2CrhZ20LO/e1di2dNRDwrbr\nQeh2fU/c9pv3NINkt9/b/Ma6+c4D4f633N3cflsSxL31rmbbKHArSd8I+rgzCPKmVW6jcG3SNq6K\nmVS6jKrUBmG8qJ2UVJ7NgqlBsDSr8BqGaysCr1GItSY017bybdrviMK1Q6icW1PRtkbbMPA8K+Jm\nbaOAbRTkzYK4a4eCirjJ+0T0vVAVzg1eQxbEPSq4SeNg8v6z07n7P+t7DFst7KQFAIAxG8BKy+Aw\naQEAYIB8jUnLVhSXAwAAo8BKCwAAA8TloaaFnbSsHLW3sW3pmGPDtmtHNLd/40AczLo1qDJ74+33\nNLZ95c57w/2/HoRr//GbcdsoXHt7EOSVpDuCce0Pth06mIVrm6G3Q0kQNwrNZgHOqG0Uxmsbgs3a\ndreiEaUAAAoQSURBVFfltjwgHO/f/s2obR9dvSEu6pK2TeJgZ/H+S93sH22vCu0mVU+jPkqr5ErS\n8mrz+2B5JX4Nhw42t0eBf6k8oLs3CeIes2elsS2qkjsETFqauDwEAABGYWFXWgAAGDNWWpqYtAAA\nMEBMWpq4PAQAAEZhYVdaouq3fvRxYds7l5qh3a/dEVeevf7WuxvbbrqjGcS96db94f5fu6MZuv16\nEtq9I6h+e+DeZrhWkg4dDB4DH4Rus3BtFMRdSx4jH4VQ21a0zUKwVRVxg79Ksr9Uav6CKQ2WVvU5\ngCDuUI+1CLoK3Za2TUO7QZi4JuAbbcsq4k7ubW7PgrjLQWg2e0+JgvzRk/xWl+NjHbm7+Wvv2L2r\nYdu+8XPXtLCTFgAAxmydSUsDl4cAAMAosNICAMAAcXmoiUkLAAADxKSlaWEnLcsPP6Gx7eBDHxW2\nvfHWZhD2szffEba97pa7Gtu+9PVmOPcr34iDuHcFFXHv3R+Haw8GodsocCvFodkoiJuGa4MgbB7E\nbRd4jYKtVft3FHitadvVteYu3qQW4Y2vpspu28q189ZFaHepZZA32x4eazkOsUb/DpPVPWHbKMx7\n4N7419PBPWWVvrPKuXuCMPC3HLM7bNu3Ra0u3QaZFgAAMAoLu9ICAMCYLcIq6awxaQEAYICYtDRx\neQgAAIzC4q60POpJjU3XBoFbSbr4S7c2tl3xpW+Ebb/4tTsb2+76ZrPfe+5uBm4l6eA9QZXaIESW\nbV87FId21w81j7d+sLltFoHXaHsWTC39S6GLCrW1/XbZxxCPVcPX41D2TmdL8/u7r21gt7aPMIjb\nsqLuJAntRmHee5O29+w5stn2iOb7YnQjgyRdE2w77uhdYdunPvKYcPu8DPX9oE+LO2kBAGDE+GOh\nictDAABgFFhpAQBggLg81MSkBQCAAWLS0rSwk5ar1x/W2PY7F38+bPv31329se2OW+OKtlHA9sBd\nzeq5awfi/aPA7FqwTSqvJpu2bV1Ntvx6atsfrq5+OMf2lFQqYCIz70q/pVV124Z7p9ubSYWs0u7S\nSnP7/qDS7p174xDt3Xc+tLHt4rCl9C+//fjkK+jLwk5aAAAYs7H90TUPTFoAABggVl6bmLQAADBA\nZFqaep20mNlZkt4oaSLpD9z9DUGb35X0/ZLulvRyd7+ipO+n/ot/P8uhAgAWwM3Btus/nDR+0W93\nOJJhMrNfkvTjkv5xY9Nr3f0D/Y3ocL3VaTGziaTzJJ0l6YmSzjazJ2xp83xJj3X3UyT9pKQ3zX2g\nO8j6HV/pewgLg3M5O5zL2eA8jo+vr83ko/awkn7L3U/b+BjMhEXqt7jcGZKuc/cvuvtBSedLeuGW\nNj8o6Q8lyd0vlbTPzB4x32HuHH7nV/sewsLgXM4O53I2OI/j09OkRZJs1q9lVvqctJwg6YZNn9+4\nse3B2pzY8bgAANjJXm1mV5rZ28xsX9+D2azPTIsXtts64yvdDwCA0eoqiGtmF0mKitC8TtMYxn/c\n+PyXJf2mpB/rZCDbYO79zAHM7OmSfsndz9r4/LWS1jeHcc3szZI+7O7nb3x+jaRnu/vNW/piIgMA\nmBt37/QSyqx/r21nvGb2aEnvc/cnzXIsbfS50nK5pFM2TsqXJb1Y0tlb2rxX0jmSzt+Y5Hxj64RF\n6v6bBwCAeerr95qZfYu735fa/iFJn+ljHJneJi3ufsjMzpF0oaa3PL/N3a82s1dtfP0t7v5+M3u+\nmV0n6S5Jr+hrvAAA7ABvMLMnaxrF+IKkV/U8nsP0dnkIAACgRp93D1Uzs7PM7Bozu9bMfj5p87sb\nX7/SzE6b9xjH4MHOo5n9yMb5+7SZfczMvrOPcY5ByffkRrvTzeyQmb1onuMbi8Kf7eeY2RVm9lkz\n+/CchzgaBT/fDzOzD5jZpzbO5ct7GObgmdnbzexmM0svj/D7pgfuPooPTS8hXSfp0ZJWJH1K0hO2\ntHm+pPdv/P/TJP1N3+Me2kfhefxfJB2z8f9ncR63fy43tfugpL+Q9K/6HvfQPgq/J/dJ+ntJJ258\n/rC+xz3Ej8Jz+UuSfu2+8yjp65KW+x770D4kPUvSaZI+k3yd3zc9fIxppYVidLPxoOfR3T/h7rdv\nfHqpqI2TKfmelKRXS3q3/qksNg5Xch5fIulP3f1GSXL3W+Y8xrEoOZdfkXT0xv8fLenr7n5ojmMc\nBXe/RNJtD9CE3zc9GNOkhWJ0s1FyHjf7MUnv73RE4/Wg59LMTtD0l8Z9j6AgRNZU8j15iqSHmtmH\nzOxyM/vRuY1uXErO5VslfbuZfVnSlZJ+ek5jWzT8vunBmJ7yTDG62Sg+H2b2vZJeKemZ3Q1n1ErO\n5Rsl/YK7u5mZBlweu0cl53FF0lMknSlpr6RPmNnfuPu1nY5sfErO5bmSPuXuzzGzkyVdZGanuvsd\nHY9tEfH7Zs7GNGm5SdJJmz4/SdOZ7QO1OXFjG/5JyXnURvj2rZLOcvcHWiLdyUrO5XdpWmdImuYH\nvt/MDrr7e+czxFEoOY83SLrF3fdL2m9mF0s6VRKTlsOVnMtnSPoVSXL3683sC5Ier2ntLJTj900P\nxnR56P5idGa2qmkxuq1v/O+V9DLp/oq7YTG6He5Bz6OZfauk90h6qbtf18MYx+JBz6W7/zN3f4y7\nP0bTXMu/ZcLSUPKz/d8lfbeZTcxsr6bBx6vmPM4xKDmX10j6PknayGA8XtLn5zrKxcDvmx6MZqXF\nKUY3EyXnUdIvSnqIpDdtrBAcdPcz+hrzUBWeSzyIwp/ta8zsA5I+LWld0lvdnUnLFoXfk78q6R1m\ndqWmf7j+nLvf2tugB8rM/ljSsyU9zMxukPR6TS9T8vumRxSXAwAAozCmy0MAAGAHY9ICAABGgUkL\nAAAYBSYtAABgFJi0AACAUWDSAgAARoFJCwAAGAUmLQAAYBSYtAAAgFEYTRl/ANtnZj+p6QMbv03S\nH0l6lKTjJD1J0zLujYdmAsDQUMYfWHBm9hOSPuXul5nZ6ZIukvRyTZ+XcqGk73f3C3scIgAUYaUF\nWHzHuvtlG///KEnr7v7nZrZH0rPd/ZIexwYAxci0AAvO3X9906fPkfSRje37t05YzOxkM3v7HIcH\nAMVYaQF2ljMlvTn6gpmdI+m7JD16ngMCgFKstAALzMwmZvZcM1sys0dKerw2Vlo2vv6a+/7f3c+T\n9M75jxIAyjBpARbbqzQN254i6cWS7pZ0oySZ2QskXbWlvc11dABQgctDwGL7mKR3aTphuVLST0n6\nDTP7oqTr3f1dPY4NAKowaQEWmLtfKelHt2z+r32MBQDa4vIQAAAYBSYtAABgFJi0AJB0f+Xc10h6\nkpn932b2uL7HBACbUcYfAACMAistAABgFJi0AACAUWDSAgAARoFJCwAAGAUmLQAAYBSYtAAAgFFg\n0gIAAEaBSQsAABgFJi0AAGAU/n9bhqlv51LtqgAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# graficando la solución\n", "u_mat = u_sol.vector().array().reshape(N1+1, N2+1)\n", "\n", "x = np.linspace(0, 1, N1+2)\n", "y = np.linspace(0, 1, N1+2)\n", "X, Y = np.meshgrid(x, y)\n", "\n", "fig, ax = plt.subplots(1, 1, figsize=(8, 6))\n", "\n", "c = ax.pcolor(X, Y, u_mat, vmin=-5, vmax=5, cmap=mpl.cm.get_cmap('RdBu_r'))\n", "cb = plt.colorbar(c, ax=ax)\n", "ax.set_xlabel(r\"$x_1$\", fontsize=18)\n", "ax.set_ylabel(r\"$x_2$\", fontsize=18)\n", "cb.set_label(r\"$u(x_1, x_2)$\", fontsize=18)\n", "fig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Para profundizar en como utilizar el *framework* [FEniCS](https://fenicsproject.org/), les recomiendo que visiten la [documentación](https://fenicsproject.org/documentation/tutorial/index.html) del sitio, que tienen varios ejemplos.\n", "\n", "Con esto concluyo este artículo. Obviamente, no es más que una introducción al fascinante y complejo mundo de las [Ecuaciones en derivadas parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales), cada clase de [EDP](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) es un mundo en sí mismo y quedaron muchos temas sin tratar; los cuales tal vez profundice en algún otro artículo. Espero que les pueda servir como referencia y lo hayan encontrado instructivo.\n", "\n", "Saludos!\n", "\n", "*Este post fue escrito utilizando Jupyter notebook. Pueden descargar este [notebook](https://github.com/relopezbriega/relopezbriega.github.io/blob/master/downloads/pyPDE.ipynb) o ver su version estática en [nbviewer](https://nbviewer.ipython.org/github/relopezbriega/relopezbriega.github.io/blob/master/downloads/pyPDE.ipynb).*" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "name": "python", "version": "2.7.11" } }, "nbformat": 4, "nbformat_minor": 0 }