Esta notebook fue creada originalmente como un blog post por Raúl E. López Briega en Matemáticas, análisis de datos y python. El contenido esta bajo la licencia BSD.
Nuestra comprensión de los procesos fundamentales de la naturaleza se basa en gran medida en Ecuaciones 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 también juegan un papel central en las matemáticas modernas, especialmente en la geometría y el análisis; lo que las convierte en una herramienta de suma utilidad que debemos conocer.
Nota: Este artículo corresponde a la tercer entrega de mi serie de artículos sobre Cálculo con Python; los anteriores fueron: Introducción al Cálculo y Ecuaciones Diferenciales con Python, los cuales es recomendable haber leído previamente.
Una Ecuación en derivadas parciales es una ecuación que, como su nombre lo indica, contiene derivadas parciales. A diferencia de lo que habíamos visto con las ecuaciones diferenciales ordinarias, en donde la función incógnita depende solo de una variable; en las Ecuaciones en derivadas parciales, o EDP 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 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 de $u$. Lo podríamos expresar de la siguiente forma:
$$F(x, y, u(x, y), u_x(x, y), u_y(x, y)) = F(x, y, u, u_x, u_y) = 0$$Al igual de como pasaba con las ecuaciones diferenciales ordinarias, el orden de una Ecuación en derivadas parciales va a estar dado por la mayor derivada presente. Por lo tanto, el caso que expresamos más arriba corresponde a una EDP, con dos variables independientes, de primer orden. Si quisiéramos expresar una EDP de segundo orden, podríamos hacerlo de la siguiente manera:
$$F(x, y, u, u_x, u_y, u_{xx}, u_{xy}, u_{yy})=0.$$La clasificación de las Ecuaciones 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:
1- El orden de la EDP: Como ya mencionamos, el mismo va a estar dado por el orden de la mayor derivada presente.
2- Número de variables: Esta clasificación va a estar dada por la cantidad de variables independientes que contenga la EDP.
3- Linealidad: Esta es una de las clasificaciones más importantes, vamos a poder clasificar a las Ecuaciones en derivadas parciales en lineales o no lineales. En las lineales, 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 de segundo orden con dos variables, va a tomar la siguiente forma:
$$Au_{xx} + Bu_{xy} + Cu_{yy} + Du_x + Eu_y + Fu = G$$en donde $A, B, C, D, E, F,$ y $G$ pueden ser constantes o una función dada de $x$ e $y$. Por ejemplo:
$u_{tt} = e^tu_{xx} + \sin t$, sería una ecuación lineal.
$uu_{xx} + u_y = 0$, sería una ecuación no lineal.
$u_{xx} + yu_{yy} + u_x = 0$, sería una ecuación lineal.
$xu_x + yu_y + u^2 = 0$, sería una ecuación no lineal.
Tipos de ecuaciones lineales: Asimismo, a las Ecuaciones en derivadas parciales lineales de segundo orden las vamos a poder subdividir en las siguientes categorías:
Ecuaciones parabólicas: Las cuales van a describir el flujo del calor y el proceso de difusión. Éstas ecuaciones van a satisfacer la condición $B^2 -4AC = 0$.
Ecuaciones hiperbólicas: Las cuales describen los sistemas de vibración y los movimientos de ondas. Satisfacen la condición $B^2 -4AC > 0$.
Ecuaciones Elípticas: Las cuales describen los fenómenos de estados estacionarios y satisfacen la condición $B^2 -4AC < 0$.
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.
5- Tipos de coeficientes: Por último, podemos clasificar a las EDP 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.
Existen varios métodos que podemos utilizar para intentar resolver las Ecuaciones 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 (EDO), o en alguna ecuación algebraica; las cuales son más sencillas de resolver. Algunos los métodos que podemos utilizar son:
1- El método de separación 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. La idea es reducir a la EDP de $n$ variables, en $n$ EDOs.
2- El método de la transformada integral: Este método es similar al que ya vimos al resolver EDOs. La idea es aplicar una transformada integral para reducir una EDP de $n$ variables, en otra de $n - 1$ variables. De esta forma una EDP de 2 variables, puede ser transformada en una EDO.
3- El método del cambio de coordenadas: Este método intenta cambiar a la EDP original en una EDO o en otra EDP más sencilla de resolver por medio del cambio de coordenadas del problema.
4- El método de perturbación: Este método aplica la teoría perturbacional para intentar cambiar un problema de EDP no lineal en una serie de problemas lineales que se aproximan al no lineal original.
5- El método de expansión de autofunciones: Este método intentan encontrar la solución de una EDP como una suma infinita de autofunciones. Estas autofunciones son halladas por medio de la resolución del problema del valor propio que corresponde al problema original.
6- El método de las ecuaciones integrales: Este método convierte a la EDP en una ecuación integral; la cual luego es resuelta aplicando ciertas técnicas particulares que se aplican a ese tipo ecuaciones.
7- Métodos numéricos: La mayoría de las técnicas para resolver numéricamente a las EDP se basan en la idea de discretizar el problema en cada variable independiente que se produce en la EDP, y de esta forma, reformular el problema en una forma algebraica. Esto usualmente resulta en problemas de álgebra lineal de gran escala. Dos de las técnicas principales para reformular las EDP a una forma algebraica son, los métodos de diferencias finitas (MDF), donde las derivadas del problema son aproximadas por medio de la fórmula de diferencias finitas; y los métodos 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.
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.
La más simple de las Ecuaciones en derivadas parciales que nos podemos encontrar es:
$u_x = 0$, en donde $u = u(x, y)$.
Esta ecuación nos dice que la 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.
Subiendo un poco más la complejidad, podemos pasar a una EDP de segundo orden, como la siguiente:
$u_{xx} = 0$
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.
Por último, si quisiéramos resolver la siguiente EDP:
$u_{xy} = 0$
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:
$u(x, y) = F(y) + g(x)$, en donde $F' = f$
Como podemos ver de estos ejemplos, las soluciones analíticas de las EDP dependen de funciones arbitrarias (en lugar de constantes arbitrarias como era el caso de las EDO). Por lo tanto vamos a necesitar condiciones auxiliares para poder determinar una única solución.
Al igual que nos pasaba cuando vimos las ecuaciones diferenciales ordinarias; las EDP 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 del problema que estemos analizando y pueden llevar a ser de dos tipos diferentes: condiciones iniciales y condiciones de frontera.
La condición inicial va a establecer el estado del problema al momento de tiempo cero, $t_0$. Por ejemplo para el problema de difusión, la condición inicial va a ser:
$$u(x, t_0) = \phi(x)$$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, $\phi(x)$ va a representar la temperatura inicial.
La condición de frontera nos va a delimitar el dominio en el que nuestra EDP es válida. Así por ejemplo, volviendo al problema de difusión, el dominio en el que nuestra EDP es válida, puede estar delimitado por la superficie del objeto que contiene al líquido. Existen varios tipos de condiciones de frontera, de las cuales las más importantes son:
La condición de frontera de Dirichlet, en dónde los valores válidos de la función incógnita $u$ son especificados.
La condición de frontera de Neumann, en donde los valores válidos especificados son dados para alguna de las derivadas de $u$.
La condición de frontera de Robin, en donde los valores válidos son especificados por una combinación lineal de una función y las derivadas de $u$.
Las EDP de primer orden poseen una interpretación geométrica la cual nos puede facilitar alcanzar una solución general para ellas.
Tomemos la siguiente ecuación de coeficientes constantes:
$au_x + bu_y = 0$, en donde $a$ y $b$ son constantes y ambas no pueden ser cero.
En esta ecuación la cantidad $au_x + bu_y$ es la derivada direccional de $u$ en la dirección del 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 $(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:
$u(x, y) = f(bx - ay)$, en donde $f$ es una función arbitraria de una variable.
Entonces, si por ejemplo, quisiéramos resolver la EDP, $4u_x - 3u_y= 0$, con la condición 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 que satisface la condición de frontera es $u(x, y) = (3x + 4y)^3 / 64$.
Consideremos ahora la siguiente EDP de coeficientes variables:
$u_x + yu_y = 0$
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:
$\frac{dy}{dx} = \frac{y}{1}$.
Podemos resolver esta ecuación como una EDO y así obtener:
$y = Ce^x$, o lo que es lo mismo $e^{-x}y = C$.
La solución de nuestra EDP entonces va a ser constante en estas curvas características (ver gráfico); y van a responder a la siguiente solución general:
$u(x, y) = f(e^{-x}y)$, en donde $f$ es una función arbitraria.
Es tiempo de nuevamente recurrir a nuestros queridos paquetes científicos de Python, NumPy, Matplotlib, SymPy y SciPy para ayudarnos a resolver las Ecuaciones en derivadas parciales. Así como en el caso de las Ecuaciones diferenciales ordinarias vimos que existía dentro del paquete SymPy, el solucionador genérico sympy.dsolve; para el caso de las EDP, vamos a tener al solucionador sympy.pdsolve. Aunque en este caso, es mucho más limitado que su versión para EDO; ya que solo vamos a poder resolver EDP de primer orden. Veamos como funciona:
In [1]:
# <!-- collapse=True -->
# importando modulos necesarios
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
import sympy
# imprimir con notación matemática.
sympy.init_printing(use_latex='mathjax')
In [2]:
# Defino las variables
x, y, u, z = sympy.symbols('x y u z')
f = sympy.Function('f')
In [3]:
# Defino la EDP 4u_x - 3u_y = 0
u = f(x, y)
u_x = u.diff(x)
u_y = u.diff(y)
eq = sympy.Eq(4*u_x - 3*u_y)
eq
Out[3]:
In [4]:
# Resuelvo la ecuación
sympy.pdsolve(eq)
Out[4]:
In [5]:
# Defino la EDP u_x + yu_y = 0
u = f(x, y)
u_x = u.diff(x)
u_y = u.diff(y)
eq2 = sympy.Eq(u_x + y*u_y)
eq2
Out[5]:
In [6]:
# Resuelvo la ecuación
sympy.pdsolve(eq2)
Out[6]:
In [7]:
# Calsificación de EDP.
sympy.classify_pde(eq2)
Out[7]:
Como podemos comprobar, obtuvimos los mismos resultados utilizando sympy.pdsolve que en nuestro análisis manual con la interpretación geométrica. Otra limitación que vamos a tener al trabajar con sympy.pdsolve 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 también nos ofrece la función classify_pde la cual nos ayuda a saber con que tipo de EDP estamos tratando. (Recordemos que pdsolve solo puede resolver EDPs de primer orden).
Otra forma en que nos podemos ayudar de SymPy para resolver EDPs, es utilizando el método de separación de variables. La característica esencial de esta técnica es transformar la EDP en un conjunto de EDOs (las cuales podemos solucionar con la ayuda de SymPy). De esta forma, la solución requerida de la EDP se expresa como 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 pueden ser resueltos por este método. Para ilustrar como funciona esta técnica, veamos un ejemplo. Vamos a resolver la siguiente EDP.
$y^2u_x^2 + x^2u_y^2 = (xyu)^2$, que cumple con la condición $u(x, 0) = 3e^{x^2/4}$.
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:
$$y^2(f'(x) g(y))^2 + x^2(f(x) g'(y))^2 = x^2 y^2 (f(x) g(y))^2$$lo que es equivalente a decir;
$$\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$$Luego, ayudándonos de la constante de separación $\lambda^2$, podemos separar a esta ecuación en dos EDOs, del siguiente modo; primero igualamos la ecuación anterior a $\lambda^2$:
$$\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$$y luego separamos ambas ecuaciones para obtener:
$$\frac{1}{x}\frac{f'(x)}{f(x)} = \lambda \\ \frac{g'(x)}{yg(y)}= \sqrt{1 - \lambda^2}$$Ahora podemos utilizar sympy.dsolve para resolver ambas EDOs:
In [8]:
# EDO n° 1
edo1 = sympy.Eq((1 / x) * (f(x).diff(x)/f(x)) - z)
edo1
Out[8]:
In [9]:
# Resolviendo EDO n° 1
sympy.dsolve(edo1)
Out[9]:
In [10]:
# EDO n° 2
edo2 = sympy.Eq((f(y).diff(y)) / (y*f(y)) - sympy.sqrt(1 - z**2))
edo2
Out[10]:
In [11]:
# Resolviendo EDO n° 2
sympy.dsolve(edo2)
Out[11]:
Entonces ahora podemos utilizar estos resultados para armar la solución final a nuestra EDP original, el cual va a ser:
$$u(x, y) = C e^{\frac{\lambda}{2}x^2 + \frac{1}{2}y^2\sqrt{1 - \lambda^2}}$$En donde $C = C_1 C_2$ es una constante arbitraria.
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:
$$u(x, y) = 3 e^{\frac{1}{4}\left(x^2 + y^2\sqrt{3}\right)}$$Si bien debemos realizar un trabajo manual previo, aun así SymPy sigue siendo de gran ayuda para facilitarnos llegar a la solución final.
Por último, para cerrar este artículo, veamos como podemos aplicar el métodos de los elementos finitos (MEF) con Python. Para esto nos vamos ayudar de la librería FEniCS, la cual es un framework para resolver numéricamente problemas generales de EDP utilizando el métodos de los elementos finitos. Para instalar esta librería en Ubuntu, pueden utilizar los siguientes comandos:
sudo add-apt-repository ppa:fenics-packages/fenics
sudo apt-get update
sudo apt-get install fenicsDeben tener en cuenta que por ahora solo funciona con Python 2. 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.
In [12]:
# importando modulos de fenics
import dolfin
import mshr
dolfin.parameters["reorder_dofs_serial"] = False
dolfin.parameters["allow_extrapolation"] = True
El problema que vamos a resolver con la ayuda de FEniCS, va a ser la siguiente EDP:
$$u_{xx} + u_{yy} = 0$$Con las siguientes condiciones de frontera:
$$u(x=0) = 3 ; \ u(x=1)=-1 ; \ u(y=0) = -5 ; \ u(y=1) = 5$$El primer paso en la solución de una EDP utilizando el métodos de los elementos finitos, es definir una malla que describa la discretización del dominio del problema. Para este caso, vamos a utilizar la función RectangleMesh que nos ofrece FEniCS.
In [13]:
# Discretizando el problema
N1 = N2 = 75
mesh = dolfin.RectangleMesh(dolfin.Point(0, 0), dolfin.Point(1, 1), N1, N2)
El siguiente paso es definir una representación del 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, 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.
In [14]:
# Funciones bases
V = dolfin.FunctionSpace(mesh, 'Lagrange', 1)
u = dolfin.TrialFunction(V)
v = dolfin.TestFunction(V)
Ahora debemos definir a nuestra EDP en su formulación débil equivalente para poder tratarla como un problema de álgebra lineal que podamos resolver con el MEF.
In [15]:
# Formulación debil de la EDP
a = dolfin.inner(dolfin.nabla_grad(u), dolfin.nabla_grad(v)) * dolfin.dx
f = dolfin.Constant(0.0)
L = f * v * dolfin.dx
Por último, solo nos falta definir las condiciones de frontera.
In [16]:
# Defino condiciones de frontera
def u0_top_boundary(x, on_boundary):
return on_boundary and abs(x[1]-1) < 1e-8
def u0_bottom_boundary(x, on_boundary):
return on_boundary and abs(x[1]) < 1e-8
def u0_left_boundary(x, on_boundary):
return on_boundary and abs(x[0]) < 1e-8
def u0_right_boundary(x, on_boundary):
return on_boundary and abs(x[0]-1) < 1e-8
In [17]:
# Definiendo condiciones de frontera de Dirichlet
bc_t = dolfin.DirichletBC(V, dolfin.Constant(5), u0_top_boundary)
bc_b = dolfin.DirichletBC(V, dolfin.Constant(-5), u0_bottom_boundary)
bc_l = dolfin.DirichletBC(V, dolfin.Constant(3), u0_left_boundary)
bc_r = dolfin.DirichletBC(V, dolfin.Constant(-1), u0_right_boundary)
# Lista de condiciones de frontera
bcs = [bc_t, bc_b, bc_r, bc_l]
Con esta especificación de las condiciones de frontera, ya estamos listos para resolver nuestra EDP utilizando la función dolfin.solve. El vector resultante, luego lo podemos convertir a una matriz de NumPy y utilizarla para graficar la solución con Matplotlib.
In [18]:
# Resolviendo la EDP
u_sol = dolfin.Function(V)
dolfin.solve(a == L, u_sol, bcs)
In [19]:
# graficando la solución
u_mat = u_sol.vector().array().reshape(N1+1, N2+1)
x = np.linspace(0, 1, N1+2)
y = np.linspace(0, 1, N1+2)
X, Y = np.meshgrid(x, y)
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
c = ax.pcolor(X, Y, u_mat, vmin=-5, vmax=5, cmap=mpl.cm.get_cmap('RdBu_r'))
cb = plt.colorbar(c, ax=ax)
ax.set_xlabel(r"$x_1$", fontsize=18)
ax.set_ylabel(r"$x_2$", fontsize=18)
cb.set_label(r"$u(x_1, x_2)$", fontsize=18)
fig.tight_layout()
Para profundizar en como utilizar el framework FEniCS, les recomiendo que visiten la documentación del sitio, que tienen varios ejemplos.
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, cada clase de EDP 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.
Saludos!
Este post fue escrito utilizando Jupyter notebook. Pueden descargar este notebook o ver su version estática en nbviewer.