Una Ecuación diferencial es una ecuación que involucra una variable dependiente y sus derivadas con respecto a una o más variables independientes. Muchas de las leyes de la naturaleza, en Física, Química, Biología, y Astronomía; encuentran la forma más natural de ser expresadas en el lenguaje de las Ecuaciones diferenciales. Estas ecuaciones no sólo tienen aplicaciones en la ciencias físicas, sino que también abundan sus aplicaciones en las ciencias aplicadas como ser Ingeniería, Finanzas y Economía.
Es fácil entender la razón detrás de esta amplia utilidad de las Ecuaciones diferenciales. Si recordamos que $y = f(x)$ es una función, entonces su derivada $dy / dx$ puede ser interpretada como el ritmo de cambio de $y$ con respecto a $x$. En muchos procesos naturales, las variables involucradas y su ritmo de cambio están conectados entre sí por medio de los principios científicos básicos que rigen el proceso. Cuando esta conexión es expresada matemáticamente, el resultado generalmente es una Ecuación diferencial.
La ley del enfriamiento de Newton o enfriamiento newtoniano establece que la tasa de pérdida de calor de un cuerpo es proporcional a la diferencia de temperatura entre el cuerpo y su ambiente.
La ley del enfriamiento de Newton hace una declaración sobre una tasa instantánea de cambio de la temperatura. Por lo que si traducimos esto al lenguaje de las matemáticas, podemos arribar a una Ecuación diferencial.
Lo que la ley nos dice, es que el ritmo de cambio de la temperatura $\frac{dT}{dt}$ es proporcional a la tempuratura de nuestro objeto $T(t)$ y la tempuratura del ambiente $T_{a}$, por lo tanto esto lo podemos expresar de la siguiente manera:
$$\frac{dT}{dt} = -k(T - T_{a})$$En donde $k$ es una constante de proporcionalidad, $t$ es la variable tiempo, $T_{a}$ es la temperutara del ambiente y $T$ es la función incognita que queremos encontrar.
A las Ecuaciones diferenciales las podemos clasificar en dos grandes grupos:
Ecuaciones diferencias ordinarias
Ecuaciones en derivadas parciales
La Ecuación diferencial de la ley del enfriamiento de Newton, es el caso típico de una Ecuación diferencial ordinaria, ya que todas las derivadas involucradas son tomadas con respecto a una única y misma variable independiente (en este caso, el tiempo).
Cuando una Ecuación diferencial contiene derivadas con respecto a única variable independiente, las denominamos Ecuaciones diferenciales ordinarias o EDO para abreviar.
Ejemplos de estas ecuaciones son:
$$\frac{d^2y}{dx^2} + \frac{dy}{dx} + y = 0; \hspace{2cm} \frac{dy}{dx} = 2xy; $$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.
Ejemplos de estas ecuaciones son:
$$\frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} + \frac{\partial^2 w}{\partial z^2} = 0 \hspace{2cm} a^2\left(\frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} + \frac{\partial^2 w}{\partial z^2}\right) = \frac{\partial w}{\partial t}$$La clasificación de las Ecuaciones diferenciales es algo muy importante, ya que dependiendo del tipo de ecuación con el que estemos tratando, distintos seran los caminos que podemos utilizar para resolverlas.
Las podemos clasificar de la siguiente manera:
El orden de una Ecuación diferencial va a ser igual al orden de la mayor derivada presente. Así, en nuestro primer ejemplo, la Ecuación diferencial de la ley del enfriamiento de Newton es de primer orden, ya que nos encontramos ante la primer derivada de la temperatura con respecto al tiempo. La ecuación general de las ecuaciones diferenciales ordinarias de grado $n$ es la siguiente:
$$F\left(x, y, \frac{dy}{dx}, \frac{d^2y}{dx^2}, \dots , \frac{d^ny}{dx^n}\right) = 0 $$o utilizando la notación prima para las derivadas,
$$F(x, y, y', y'', \dots, y^{(n)}) = 0$$La más simple de todas las ecuaciones diferenciales ordinarias es la siguiente ecuación de primer orden:
$$ \frac{dy}{dx} = f(x)$$y para resolverla simplemente debemos calcular su integral indefinida:
$$y = \int f(x) dx + c$$.
Una ecuación separable es una ecuación diferencial de primer orden en la que la expresión para $dx / dy$ se puede factorizar como una función de x multiplicada por una función de y. En otras palabras, puede ser escrita en la forma:
$$\frac{dy}{dx} = f(x)g(y)$$El nombre separable viene del hecho de que la expresión en el lado derecho se puede "separar" en una función de $x$ y una función de $y$.
Para resolver este tipo de ecuaciones, podemos reescribirlas en la forma diferencial:
$$\frac{dy}{g(y)} = f(x)dx$$y luego podemos resolver la ecuación original integrando:
$$\int \frac{dy}{g(y)} = \int f(x) dx + c$$Éstas suelen ser las Ecuaciones diferenciales más fáciles de resolver, ya que el problema de resolverlas puede ser reducido a un problema de integración; a pesar de que igualmente muchas veces estas integrales pueden ser difíciles de calcular.
Uno de los tipos más importantes de Ecuaciones diferenciales son las Ecuaciones diferenciales lineales. Este tipo de ecuaciones son muy comunes en varias ciencias y tienen la ventaja de que pueden llegar a ser resueltas en forma analítica ya que su ecuación diferencial de primer orden adopta la forma:
$$\frac{dy}{dx} + P(x)y = Q(x) $$donde, $P$ y $Q$ son funciones continuas de $x$ en un determinado intervalo. Para resolver este tipo de ecuaciones de la forma $ y' + P(x)y = Q(x)$, debemos multiplicar los dos lados de la ecuación por el factor de integración $e^{\int P(x) dx}$ y luego integrar ambos lados.
Otras clasificaciones que comúnmente se utilizan son:
Esta clasificación va a estar dada por la cantidad de variables independientes que contenga la Ecuación diferencial
Aquí clasificamos a la Ecuación diferencial según sus coeficientes, si los mismos son constantes, se dice que la ecuación es de coeficientes constantes, en caso contrario será de coeficientes variables.
Una Ecuación diferencial homogénea es una ecuación igual a cero en la que solo hay derivadas de $y$ y términos $y$, como por ejemplo la siguiente ecuación:
$$\frac{d^4y}{dx^4} + \frac{d^2y}{dx^2} + y^2 = 0$$Las Ecuaciones diferenciales 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$.
Cuando comenzamos a lidiar con las Ecuaciones diferenciales, veremos que existen un gran número de ellas que no pueden ser resueltas en forma analítica utilizando los principios del Cálculo integral y el Cálculo diferencial; pero sin embargo, tal vez podamos encontrar soluciones aproximadas para este tipo de ecuaciones en términos de Series de potencias.
Una Serie de potencias es una serie, generalmente infinita, que posee la siguiente forma:
$$\sum_{n=0}^{\infty} C_nX^n = C_0 + C_1X + C_2X^2 + C_3X^3 + \dots $$En dónde $X$ es una variable y las $C_n$ son constantes o los coeficientes de la serie. Una Serie de potencias puede converger para algunos valores de $X$ y divergir para otros valores de $X$. La suma de la serie es una función.
$$f(x) = C_0 + C_1X + C_2X^2 + C_3X^3 + \dots + C_nX^n + \dots$$El dominio de esta función va a estar dado por el conjunto de todos los $X$ para los que la serie converge.
Las Series de Taylor son un caso especial de Series de potencias cuyos términos adoptan la forma $(x - a)^n$. Las Series de Taylor nos van a permitir aproximar funciones continuas que no pueden resolverse en forma analítica y se van a calcular a partir de las derivadas de estas funciones. Su definición matemática es la siguiente:
$$f(x) = \sum_{n=0}^{\infty}\frac{f^{(n)}(a)}{n!}(x - a)^n$$Lo que es equivalente a decir:
$$f(x) = f(a) + \frac{f'(a)}{1!}(x -a) + \frac{f''(a)}{2!}(x -a)^2 + \frac{f'''(a)}{3!}(x -a)^3 + \dots$$Una de las razones de que las Series de Taylor sean importantes es que nos permiten integrar funciones que de otra forma no podíamos manejar.
Las Series de Fourier son series infinitas expresadas en términos de seno y coseno que convergen en una función periódica y continua. Así por ejemplo una función con un período de $2\pi$, va a adaptar la forma de la siguiente serie:
$$f(x) = \frac{a_{0}}{2} + \sum_{k=1}^{\infty}(a_{k} \cos kx + b_{k} \sin kx)$$El análisis de Fourier suele ser una herramienta muy útil para resolver Ecuaciones diferenciales.
SymPy es una librería de Python para matemática simbólica. Tiene como objetivo convertirse en un sistema completo de álgebra computacional.
Algunas de las características que posee son:
Supongamos que hubo un asesinato y la policía lleva a la escena del crimen a las 2:41 am. El forense toma la temperatura de la víctima y encuentra que es de 34.5° C; luego vuelve a tomar la temperatura una hora después y la lectura marca 33.9° C.
Suponiendo que ese día el clima estuvo templado y el ambiente hubo una temperatura constante de 15° C. ¿Cuál fue la hora del asesinato?
In [1]:
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
import sympy
from scipy import integrate
# imprimir con notación matemática.
sympy.init_printing(use_latex='mathjax')
# importando modulos de fenics
import dolfin
import mshr
%matplotlib inline
dolfin.parameters["reorder_dofs_serial"] = False
dolfin.parameters["allow_extrapolation"] = True
In [2]:
def plot_direction_field(x, y_x, f_xy, x_lim=(-5, 5), y_lim=(-5, 5), ax=None):
"""Esta función dibuja el campo de dirección de una EDO"""
f_np = sympy.lambdify((x, y_x), f_xy, modules='numpy')
x_vec = np.linspace(x_lim[0], x_lim[1], 20)
y_vec = np.linspace(y_lim[0], y_lim[1], 20)
if ax is None:
_, ax = plt.subplots(figsize=(4, 4))
dx = x_vec[1] - x_vec[0]
dy = y_vec[1] - y_vec[0]
for m, xx in enumerate(x_vec):
for n, yy in enumerate(y_vec):
Dy = f_np(xx, yy) * dx
Dx = 0.8 * dx**2 / np.sqrt(dx**2 + Dy**2)
Dy = 0.8 * Dy*dy / np.sqrt(dx**2 + Dy**2)
ax.plot([xx - Dx/2, xx + Dx/2],
[yy - Dy/2, yy + Dy/2], 'b', lw=0.5)
ax.axis('tight')
ax.set_title(r"$%s$" %
(sympy.latex(sympy.Eq(y(x).diff(x), f_xy))),
fontsize=18)
return ax
In [3]:
def laplace_transform_derivatives(e):
"""
Evalua las transformadas de Laplace de derivadas de funciones sin evaluar.
"""
if isinstance(e, sympy.LaplaceTransform):
if isinstance(e.args[0], sympy.Derivative):
d, t, s = e.args
n = len(d.args) - 1
return ((s**n) * sympy.LaplaceTransform(d.args[0], t, s) -
sum([s**(n-i) * sympy.diff(d.args[0], t, i-1).subs(t, 0)
for i in range(1, n+1)]))
if isinstance(e, (sympy.Add, sympy.Mul)):
t = type(e)
return t(*[laplace_transform_derivatives(arg) for arg in e.args])
return e
Para resolver el problema, debemos utilizar la La Ecuación diferencial de la ley del enfriamiento de Newton. Los datos que tenemos son:
In [4]:
# defino las incognitas
t, k = sympy.symbols('t k')
y = sympy.Function('y')
# expreso la ecuacion
f = k*(y(t) -15)
sympy.Eq(y(t).diff(t), f)
Out[4]:
In [5]:
# Resolviendo la ecuación
edo_sol = sympy.dsolve(y(t).diff(t) - f)
edo_sol
Out[5]:
Ahora que tenemos la solución de la Ecuación diferencial, despejemos constante de integración utilizando la condición inicial.
In [6]:
# Condición inicial
ics = {y(0): 34.5}
C_eq = sympy.Eq(edo_sol.lhs.subs(t, 0).subs(ics), edo_sol.rhs.subs(t, 0))
C_eq
Out[6]:
In [7]:
C = sympy.solve(C_eq)[0]
C
Out[7]:
Ahora que ya sabemos el valor de C, podemos determinar el valor de $k$.
In [8]:
eq = sympy.Eq(y(t), C * sympy.E**(k*t) +15)
eq
Out[8]:
In [9]:
ics = {y(1): 33.9}
k_eq = sympy.Eq(eq.lhs.subs(t, 1).subs(ics), eq.rhs.subs(t, 1))
kn = round(sympy.solve(k_eq)[0], 4)
kn
Out[9]:
Ahora que ya tenemos todos los datos, podemos determinar la hora aproximada de la muerte.
In [10]:
hmuerte = sympy.Eq(37, 19.5 * sympy.E**(kn*t) + 15)
hmuerte
Out[10]:
In [11]:
t = round(sympy.solve(hmuerte)[0],2)
t
Out[11]:
In [12]:
h, m = divmod(t*-60, 60)
print "%d horas, %d minutos" % (h, m)
Es decir, que pasaron aproximadamente 3 horas y 51 minutos desde que ocurrió el crimen, por lo tanto la hora del asesinato debio haber sido alredor de las 10:50 pm.
Un método alternativo que podemos utilizar para resolver en forma analítica Ecuaciones diferenciales ordinarias complejas, es utilizar la Transformada de Laplace, que es un tipo particular de transformada integral. La idea es que podemos utilizar esta técnica para transformar nuestra Ecuación diferencial en algo más simple, resolver esta ecuación más simple y, a continuación, invertir la transformación para recuperar la solución a la Ecuación diferencial original.
Para poder comprender la Transformada de Laplace, primero debemos revisar la definición general de la transformada integral, la cuál adapta la siguiente forma:
$$T(f(t)) = \int_{\alpha}^{\beta} K (s, t) \ f(t) \ dt = F(s) $$En este caso, $f(t)$ es la función que queremos transformar, y $F(s)$ es la función transformada. Los límites de la integración, $\alpha$ y $\beta$, pueden ser cualquier valor entre $-\infty$ y $+\infty$ y $K(s, t)$ es lo que se conoce como el núcleo o kernel de la transformada, y podemos elegir el kernel que nos plazca. La idea es poder elegir un kernel que nos dé la oportunidad de simplificar la Ecuación diferencial con mayor facilidad.
Si nos restringimos a Ecuaciones diferenciales con coeficientes constantes, entonces un kernel que resulta realmente útil es $e^{-st}$, ya que al diferenciar este kernel con respecto de $t$, terminamos obteniendo potencias de $s$, que podemos equiparar a los coeficientes constantes. De esta forma, podemos arribar a la definición de la Transformada de Laplace:
$$\mathcal{L}\{f(t)\}=\int_0^{\infty} e^{-st} \ f(t) \ dt$$La principal ventaja de utilizar Transformadas de Laplace es que cambia la Ecuación diferencial en una ecuación algebraica, lo que simplifica el proceso para calcular su solución. La única parte complicada es encontrar las transformaciones y las inversas de las transformaciones de los varios términos de la Ecuación diferencial que queramos resolver. Aquí es donde nos podemos ayudar de SymPy.
Vamos a intentar resolver la siguiente ecuación:
$$y'' + 3y' + 2y = 0$$con las siguientes condiciones iniciales: $y(0) = 2$ y $y'(0) = -3$
In [13]:
# Ejemplo de transformada de Laplace
# Defino las incognitas
t = sympy.symbols("t", positive=True)
y = sympy.Function("y")
# simbolos adicionales.
s, Y = sympy.symbols("s, Y", real=True)
# Defino la ecuación
edo = y(t).diff(t, t) + 3*y(t).diff(t) + 2*y(t)
sympy.Eq(edo)
Out[13]:
In [14]:
# Calculo la transformada de Laplace
L_edo = sympy.laplace_transform(edo, t, s, noconds=True)
L_edo_2 = laplace_transform_derivatives(L_edo)
# reemplazamos la transfomada de Laplace de y(t) por la incognita Y
# para facilitar la lectura de la ecuación.
L_edo_3 = L_edo_2.subs(sympy.laplace_transform(y(t), t, s), Y)
sympy.Eq(L_edo_3)
Out[14]:
Aquí ya logramos convertir a la Ecuación diferencial en una ecuación algebraica. Ahora podemos aplicarle las condiciones iniciales para resolverla.
In [15]:
# Definimos las condiciones iniciales
ics = {y(0): 2, y(t).diff(t).subs(t, 0): -3}
ics
Out[15]:
In [16]:
# Aplicamos las condiciones iniciales
L_edo_4 = L_edo_3.subs(ics)
# Resolvemos la ecuación y arribamos a la Transformada de Laplace
# que es equivalente a nuestra ecuación diferencial
Y_sol = sympy.solve(L_edo_4, Y)
Y_sol
Out[16]:
In [17]:
# Por último, calculamos al inversa de la Transformada de Laplace que
# obtuvimos arriba, para obtener la solución de nuestra ecuación diferencial.
y_sol = sympy.inverse_laplace_transform(Y_sol[0], s, t)
y_sol
Out[17]:
In [18]:
# Comprobamos la solución.
y_sol.subs(t, 0), sympy.diff(y_sol).subs(t, 0)
Out[18]:
Las Transformadas de Laplace, pueden ser una buena alternativa para resolver Ecuaciones diferenciales en forma analítica. Pero aún así, siguen existiendo ecuaciones que se resisten a ser resueltas por medios analíticos, para estos casos, debemos recurrir a los métodos numéricos.
Supongamos ahora que queremos resolver con SymPy la siguiente Ecuación diferencial:
$$\frac{dy}{dx} = x^2 + y^2 -1$$con una condición inicial de $y(0) = 0$.
Si aplicamos lo que vimos hasta ahora, vamos a obtener el siguiente resultado:
In [19]:
# Defino incognitas
x = sympy.symbols('x')
y = sympy.Function('y')
# Defino la función
f = y(x)**2 + x**2 -1
# Condición inicial
ics = {y(0): 0}
# Resolviendo la ecuación diferencial
edo_sol = sympy.dsolve(y(x).diff(x) - f, ics=ics)
edo_sol
Out[19]:
El resultado que nos da SymPy, es una aproximación con Series de potencias (una serie de Taylor); y el problema con las Series de potencias es que sus resultados sólo suelen válidos para un rango determinado de valores. Una herramienta que nos puede ayudar a visualizar el rango de validez de una aproximación con Series de potencias son los Campos de direcciones.
Los Campos de direcciones es una técnica sencilla pero útil para visualizar posibles soluciones a las ecuaciones diferenciales de primer orden. Se compone de líneas cortas que muestran la pendiente de la función incógnita en el plano x-y. Este gráfico se puede producir fácilmente debido a que la pendiente de $y(x)$ en los puntos arbitrarios del plano x-y está dada por la definición misma de la Ecuación diferencial ordinaria:
$$\frac{dy}{dx} = f(x, y(x))$$Es decir, que sólo tenemos que iterar sobre los valores $x$ e $y$ en la grilla de coordenadas de interés y evaluar $f(x, y(x))$ para saber la pendiente de $y(x)$ en ese punto. Cuantos más segmentos de líneas trazamos en un Campo de dirección, más clara será la imagen. La razón por la cual el gráfico de Campos de direcciones es útil, es que las curvas suaves y continuos que son tangentes a las líneas de pendiente en cada punto del gráfico, son las posibles soluciones a la Ecuación diferencial ordinaria.
Por ejemplo, el Campos de direcciones de la ecuación:
$$\frac{dy}{dx} = x^2 + y^2 -1$$es el siguiente:
In [20]:
# grafico de campo de dirección
fig, axes = plt.subplots(1, 1, figsize=(7, 5))
campo_dir = plot_direction_field(x, y(x), f, ax=axes)
Ahora que ya conocemos a los Campos de direcciones, volvamos a la solución aproximada con Series de potencias que habiamos obtenido anteriormente. Podemos graficar esa solución en el Campos de direcciones, y compararla con una solución por método númericos.
En el panel de la izquierda podemos ver el gráfico de la solución aproximada por la Serie de potencias. La solución aproximada se alinea bien con el campo de direcciones para los valores de $x$ entre $-1.5$ y $1.5$, luego comienza a desviarse, lo que nos indica que la solución aproximada ya no sería válida.
In [21]:
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
# panel izquierdo - solución aproximada por Serie de potencias
plot_direction_field(x, y(x), f, ax=axes[0])
x_vec = np.linspace(-3, 3, 100)
axes[0].plot(x_vec, sympy.lambdify(x, edo_sol.rhs.removeO())(x_vec),
'b', lw=2)
# panel derecho - Solución por método iterativo
plot_direction_field(x, y(x), f, ax=axes[1])
x_vec = np.linspace(-1, 1, 100)
axes[1].plot(x_vec, sympy.lambdify(x, edo_sol.rhs.removeO())(x_vec),
'b', lw=2)
# Resolviendo la EDO en forma iterativa
edo_sol_m = edo_sol_p = edo_sol
dx = 0.125
# x positivos
for x0 in np.arange(1, 2., dx):
x_vec = np.linspace(x0, x0 + dx, 100)
ics = {y(x0): edo_sol_p.rhs.removeO().subs(x, x0)}
edo_sol_p = sympy.dsolve(y(x).diff(x) - f, ics=ics, n=6)
axes[1].plot(x_vec, sympy.lambdify(x, edo_sol_p.rhs.removeO())(x_vec),
'r', lw=2)
# x negativos
for x0 in np.arange(1, 5, dx):
x_vec = np.linspace(-x0-dx, -x0, 100)
ics = {y(-x0): edo_sol_m.rhs.removeO().subs(x, -x0)}
edo_sol_m = sympy.dsolve(y(x).diff(x) - f, ics=ics, n=6)
axes[1].plot(x_vec, sympy.lambdify(x, edo_sol_m.rhs.removeO())(x_vec),
'r', lw=2)
SciPy es un conjunto de paquetes donde cada uno ellos ataca un problema distinto dentro de la computación científica y el análisis numérico. Algunos de los paquetes que incluye, son:
scipy.integrate: que proporciona diferentes funciones para resolver problemas de integración numérica.scipy.linalg: que proporciona funciones para resolver problemas de álgebra lineal.scipy.optimize: para los problemas de optimización y minimización.scipy.signal: para el análisis y procesamiento de señales.scipy.sparse: para matrices dispersas y solucionar sistemas lineales dispersosscipy.stats: para el análisis de estadística y probabilidades.Para resolver las Ecuaciones diferenciales, el paquete que nos interesa es scipy.integrate.
SciPy nos ofrece dos solucionadores de ecuaciones diferenciales ordinarias, integrate.odeint y integrate.ode. La principal diferencia entre ambos, es que integrate.ode es más flexible, ya que nos ofrece la posibilidad de elegir entre distintos solucionadores; aunque integrate.odeint es más fácil de utilizar.
Tratemos de resolver la siguiente ecuación:
$$\frac{dy}{dx} = x + y^2$$
In [22]:
# Defino la función
f = y(x)**2 + x
f
Out[22]:
In [23]:
# la convierto en una función ejecutable
f_np = sympy.lambdify((y(x), x), f)
# Definimos los valores de la condición inicial y el rango de x sobre los
# que vamos a iterar para calcular y(x)
y0 = 0
xp = np.linspace(0, 1.9, 100)
# Calculando la solución numerica para los valores de y0 y xp
yp = integrate.odeint(f_np, y0, xp)
# Aplicamos el mismo procedimiento para valores de x negativos
xn = np.linspace(0, -5, 100)
yn = integrate.odeint(f_np, y0, xn)
Los resultados son dos matrices unidimensionales de NumPy $yp$ y $yn$, de la misma longitud que las correspondientes matrices de coordenadas $xp$ y $xn$, que contienen las soluciones numéricas de la ecuación diferencial ordinaria para esos puntos específicos. Para visualizar la solución, podemos graficar las matrices $yp$ y $yn$, junto con su Campo de direcciones.
In [24]:
# graficando la solucion con el campo de direcciones
fig, axes = plt.subplots(1, 1, figsize=(8, 6))
plot_direction_field(x, y(x), f, ax=axes)
axes.plot(xn, yn, 'b', lw=2)
axes.plot(xp, yp, 'r', lw=2)
plt.show()
En este ejemplo, solucionamos solo una ecuación. Generalmente, la mayoría de los problemas se presentan en la forma de sistemas de ecuaciones diferenciales ordinarias, es decir, que incluyen varias ecuaciones a resolver. Para ver como podemos utilizar a integrate.odeint para resolver este tipo de problemas, consideremos el siguiente sistema de ecuaciones diferenciales ordinarias, conocido el atractor de Lorenz:
Estas ecuaciones son conocidas por sus soluciones caóticas, que dependen sensiblemente de los valores de los parámetros $\sigma$, $\rho$ y $\beta$. Veamos como podemos resolverlas con la ayuda de Python.
In [25]:
# Definimos el sistema de ecuaciones
def f(xyz, t, sigma, rho, beta):
x, y, z = xyz
return [sigma * (y - x),
x * (rho - z) - y,
x * y - beta * z]
# Asignamos valores a los parámetros
sigma, rho, beta = 8, 28, 8/3.0
# Condición inicial y valores de t sobre los que calcular
xyz0 = [1.0, 1.0, 1.0]
t = np.linspace(0, 25, 10000)
# Resolvemos las ecuaciones
xyz1 = integrate.odeint(f, xyz0, t, args=(sigma, rho, beta))
xyz2 = integrate.odeint(f, xyz0, t, args=(sigma, rho, 0.6*beta))
xyz3 = integrate.odeint(f, xyz0, t, args=(2*sigma, rho, 0.6*beta))
In [26]:
# Graficamos las soluciones
from mpl_toolkits.mplot3d.axes3d import Axes3D
fig, (ax1,ax2,ax3) = plt.subplots(1, 3, figsize=(12, 4),
subplot_kw={'projection':'3d'})
for ax, xyz, c in [(ax1, xyz1, 'r'), (ax2, xyz2, 'b'), (ax3, xyz3, 'g')]:
ax.plot(xyz[:,0], xyz[:,1], xyz[:,2], c, alpha=0.5)
ax.set_xlabel('$x$', fontsize=16)
ax.set_ylabel('$y$', fontsize=16)
ax.set_zlabel('$z$', fontsize=16)
ax.set_xticks([-15, 0, 15])
ax.set_yticks([-20, 0, 20])
ax.set_zticks([0, 20, 40])
Los casos que vimos hasta ahora, se trataron de ecuaciones diferenciales ordinarias, pero ¿cómo podemos hacer para resolver ecuaciones en derivadas parciales?
Estas ecuaciones son mucho más difíciles de resolver, pero podes recurrir a la poderosa herramienta que nos proporciona el Método de los Elementos Finitos para resolverlas en forma numérica.
La idea general detrás del Método de los Elementos Finitos es la división de un continuo en un conjunto de pequeños elementos interconectados por una serie de puntos llamados nodos. Las ecuaciones que rigen el comportamiento del continuo regirán también el del elemento. De esta forma se consigue pasar de un sistema continuo (infinitos grados de libertad), que es regido por una ecuación diferencial o un sistema de ecuaciones diferenciales, a un sistema con un número de grados de libertad finito cuyo comportamiento se modela por un sistema de ecuaciones, lineales o no.
Por ejemplo, en siguiente imagen, podemos ver que en primer lugar tenemos una placa con un hueco en el centro, supongamos que queremos determinar su distribución de temperatura. Para realizar esto, deberíamos resolver la ecuación del calor para cada punto en la placa. El enfoque que utiliza el Método de los Elementos Finitos es el de dividir al objeto en elementos finitos conectados entre sí por nodos; como lo muestran la tercera y cuarta imagen. Este nuevo objeto, constituido por los elementos finitos (los triángulos de la segunda imagen) se llama malla y es una representación aproximada del objeto original. Mientras más nodos tengamos, más exacta será la solución.
El proyecto FEniCS es un framework para resolver numéricamente problemas generales de ecuaciones en derivadas parciales utilizando el métodos de los elementos finitos.
Podemos instalarlo en Ubuntu con los siguientes comandos:
sudo add-apt-repository ppa:fenics-packages/fenics
sudo apt-get update
sudo apt-get install fenicsLa 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. Por ahora solo funciona con Python 2.
El problema que vamos a resolver con la ayuda de FEniCS, va a ser la ecuación del calor en dos dimensiones en estado estacionario, definida por:
$$u_{xx} + u_{yy} = f$$donde f es la función fuente y donde tenemos 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 [2]:
# Discretizando el problema
N1 = N2 = 75
mesh = dolfin.RectangleMesh(dolfin.Point(0, 0), dolfin.Point(1, 1), N1, N2)
# grafico de la malla.
dolfin.RectangleMesh(dolfin.Point(0, 0), dolfin.Point(1, 1), 10, 10)
Out[2]:
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 [3]:
# 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 [4]:
# 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
Y definimos las condiciones de frontera.
In [5]:
# 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 [6]:
# 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]
Ahora ya podemos resolver la 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 [7]:
# Resolviendo la EDP
u_sol = dolfin.Function(V)
dolfin.solve(a == L, u_sol, bcs)
# 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)
In [8]:
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()
https://relopezbriega.github.io/
Twitter: @relopezbriega
Slides: https://relopezbriega.github.io/ecuaciones-diferenciales.html