Una de las tareas que más tiempo lleva a la hora de realizar una investigación es lo que se conoce como estudios paramétricos, para observar como influyen los diferentes valores que gobiernan el sistema analizado. Si se trabaja con hojas de cálculo, y se quiere observar como afecta cierto parámetro al sistema, tenemos que localizar la celda correspondiente e ir cambiando manualmente el valor hasta observar lo que sucede. Esto puede resultar muy tedioso, sobretodo cuando se necesita observar rápidamente como influyen varios parámetros en la evolución del sistema. IPython Notebook ofrece la posiblidad de añadir Widgets interactivos, que permiten variar de manera muy sencilla los parámetros que queramos en la manera que necesitemos.
Vamos a verlo con unos ejemplos.
En 1963 el matemático y meteorólogo del MIT, Edward Lorenz, derivó el primer sistema caótico que governaba la atmósfera. Una buena simplificación de la colección de ecuaciones fluido dinámicas proporcionan el siguiente modelo matemático;
$$ \begin{aligned} \dot{x} & = \sigma(y-x) \\ \dot{y} & = \rho x - y - xz \\ \dot{z} & = -\beta z + xy \end{aligned} $$Este es uno de los sistemas clásicos de ecuaciones diferenciales no-lineales. Exibe un gama de diferentes comportamientos para los diferentes parámetros: $\sigma$ (número de Prandt), $\beta$ (factor geométrico), $\rho$ (número de Rayleig)
Los valores normales para estos parámetros son: $\sigma$ = 10, $\beta$ = 2/3, y $\rho$ es variable, pero para un valor de 28 el sistema muestra un comportamiento caótico.
La solución del sistema para el estado estacionario es:
$$ y^s = (\rho - 1, n, n) $$Dónde:
$$ n = {\sqrt{\beta (\rho-1)}} $$Este punto de es muy inestable para $\rho$ = 28 por lo que cualqueir pequeña perturbación en las condiciones iniciales puede provocar grandes alteraciones en el sistema.
Vamos a resolver el sistema de ecuaciones anterior y a observar, con un widget interactivo, como afectan las perturbacions en las condiciones iniciales. Para ello, partiremos de la solución en estado estacionario, añadiendo la perturbación deseada.
$$ y_0 = (\rho - 1, n, n + p) $$Dónde p es la perturbación introducida al sistema, que variaremos desde 0 hasta 1.
¡¡Ok, vamos a ello!!
In [1]:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from scipy import integrate
from IPython.html.widgets import interact, interactive, fixed
from IPython.html import widgets
En orden, hemos importado:
Recuerda que en Python necesitamos cargar las librerías necesarias antes de usarlas. Otros lenguajes como MATLAB hacen esto nada más iniciarse y de forma oculta al usuario.
In [2]:
def f(y,t, beta, sigma, rho):
x1 = -beta*y[0] + y[1]*y[2]
x2 = sigma*(y[2]-y[1])
x3 = (rho - y[0])*y[1]-y[2]
return [x1, x2, x3]
Vamos a resolver de manera matricial, por lo que del modelo matemático se entiende:
$$ \\z = y[0] \\ \\x = y[1 ]\\ \\y = y[2] \\ $$Es importante recordar, que en Python comenzamos la numeración en 0.
Para poder implementar los widgets, es necesario que tengamos definida una función con los parámetros que más tarde
variaremos.
En nuestro caso, la función se llamará comportamiento_caotico, y el parámetro que cambiaremos será la perturbación introducida
al sistema.
En el apartado anteiro tan solo se ha definido el modelo matemático en una función f, pero no se ha resuelto nada, como se puede ver, no se ha introducido ningún valor. Vamos ahora a resolver el sistema en el interior de otra función. Esta nueva función tiene como parámetro la perturbación, por lo que antes de resolver el conjunto tendremos que indicar el valor de este parámetro, inicialmente tendrá un valor de 0.
In [3]:
def comportamiento_caotico(p = 0.0):
t = np.linspace(0.,100,1000)
# Parámetros del sistema en el punto definido anteriormente. Es importante definir estos parámetros, pues son
# los requeridos por la función f del punto anterior
beta = 8.0/3.0
sigma = 10.0
rho = 28.0
# Definimos n como se ha descrito arriba
n = np.sqrt(beta*(rho-1))
# Pasamos a definir las condiciones iniciales del sistema
y0 = np.array([rho-1, n, n + p])
# Y ahora mandamos resolver el sistema de ecuaciones
x = integrate.odeint(f,y0,t, args = (beta, sigma, rho) )
# Para poder hacer la representación de los datos tenemos que recuperar la información calculada con el comando anterior
# para ello es suficiente como llamar a la variables definidas en f y asignarles el valor de la matriz de resultados x.
y1 = x[:,0]
y2 = x[:,1]
y3 = x[:,2]
# Para obtener la figura, primero definimos sus dimensiones, en este caso 20x5, y solicitamos que represente en el mismo
# eje de tiempo los resultados obtenidos.
fig = plt.figure(figsize=(20, 5))
plt.plot(t,y1,'r')
plt.plot(t,y2,'b')
plt.plot(t,y3,'g')
plt.xlabel('tiempo')
plt.ylabel(u'variación(y)')
plt.show()
Fíjate que no ha pasado nada hasta ahora ya que lo único que hemos hecho ha sido crear dos funciones.
El tipo de widget que crearemos con el comando interact consite en un botón en una barra horizontal, que podremos ir moviendo de izquierda a derecha para modificar el valor de la perturbación. El comando necesita una función a la que aplicar el widget, además de conocer el parámetro a variar (o parámetros, siempre y cuando estén definidos como tal en la definición de la función), y en qué medida a de hacerlo. En nuestro caso, la función que introduciremos será la anteriro, efecto_mariposa, y como parámetro a variar la perturbación, comprendida entre 0 y 1. Recordemos la condición inicial:
$$ y_0 = (\rho - 1, n, n + p) $$
In [4]:
interact(comportamiento_caotico, p = (0.0, 1, 0.01))
Como es de esperar en un sistema caótico, pequeñas variaciones en las condiciones iniciales provocan grandes cambios en el sistema.
Pero este sistema de ecuaciones no es famoso por lo mostrado arriba. Si se realiza una representación 3D de los resultados se observan una serie de trayectorias encerradas en lo que se conoce como atractor de Lorenz.
Un atractor es el conjunto al que el sistema evoluciona después de un tiempo suficientemente largo. Para que el conjunto sea un atractor, las trayectorias que le sean suficientemente próximas han de permanecer próximas incluso si son ligeramente perturbadas. Geométricamente, un atractor puede ser un punto, una curva, una variedad o incluso un conjunto complicado de estructura fractal conocido como atractor extraño. La descripción de atractores de sistemas dinámicos caóticos ha sido uno de los grandes logros de la teoría del caos. La trayectoria del sistema dinámico en el atractor no tiene que satisfacer ninguna propiedad especial excepto la de permanecer en el atractor; puede ser periódica, caótica o de cualquier otro tipo.
Realmente, el sistema de ecuaciones fluido dinámicas mostradas al comienzo conforman en si el atractor de Lorenz. Estas ecuaciones, definen las trayectorias en el interior de un espacio de fuerzas, y para los valores estudiados en el apartado anterior se observan unas trayectorias que bien podrían tener forma de alas de mariposa, de ahí el nombre: efecto mariposa.
Para diferentes valores a los anteriormente empleados, las formas que proporciona se ven alteradas, llegando incluso a disponer las trayectorias en una forma tórica, concretamente para valores de $\beta$ cercanos a 100.
Vamos a observar primero como quedaría representado el sistema anterior en 2 dimensiones.
Primero intorduciremos un par de bibliotecas adicionales, que nos permitirán representar las gráficas de una manera vistosa.
In [5]:
from IPython.html.widgets import interact, interactive
from IPython.display import clear_output, display, HTML
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import cnames
from matplotlib import animation
La resolución del sistema de equaciones en este caso es algo más complicada, por este motivo vamos a cambiar el código empleado en el apartado anterior. La filosofía de resolución es la misma, pero con pequeños matices. Vamos a resolver de manera simultánea el sistema para diferentes valores iniciales. Vamos a trabar con todo el código de una vez, es decir, el modelo matemático y su resolución estarán directamente en el interior de la función que llamará el comando interact. Veamos cómo:
In [6]:
# En primer lugar definimos la función que llamará el comando interact.
# Tenemos que decirle los parámetros que queremos que modifique.
# N = número de valores iniciales que queremos, es decir, las veces que resolvemos simultáneamente el sistema para diferentes
# valores iniciales.
# angle = es el ángulo con el que se proyectará la figura, podremos modificarlos para observar la figura en 360º horizontales.
# max_time = es el tiempo máximo de iteración.
# sigma, beta y rho son los parámetros del sistema.
def efecto_mariposa(N=10, angle=0.0, max_time=1.0, sigma=10.0, beta=8.0/3.0 , rho=28.0):
fig = plt.figure() #definimos la figura
ax = fig.add_axes([0, 0, 1, 1], projection='3d') # definimos los ejes e indicamos que queremos una proyección 3D
ax.axis('off')
# Preparamos los límites de los ejes, prueba a modificarlos más tarde.
ax.set_xlim((-25, 25))
ax.set_ylim((-35, 35))
ax.set_zlim((5, 55))
# Definimos de nuevo el modelo matemático. Si te fijas bien observarás que es el mismo que el empleado para el apartado
# anterior, solo que ahora trabajamos directamente con las variables x, y, z. Esto lo hacemos así para no dificultar la
# tarea. En la resolución del modelo trabajaremos de manera vectorial, si empleamos aquí la misma forma que para el
# apartado anterior podemos tener problemas de convergencia, además de complicar mucho el código, al tener que trabajar
# con "matrices sobre matrices".
def lorenz_deriv((x, y, z), t0, sigma=sigma, beta=beta, rho=rho):
"""Calcula la derivada respecto del tiempo para el sistema de Lorenz"""
#Esta definición de arriba entre comillas sirve para documentar la función
dxdt = sigma * (y - x)
dydt = x * (rho - z) - y
dzdt = x * y - beta * z
return [dxdt, dydt, dzdt]
# Para definir los valores iniciales lo haremos de manera aleatoria, uniformemente distribuidos entre -15 y 15,
# para ello empleamos el comando np.random
np.random.seed(1)
x0 = -15 + 30 * np.random.random((N, 3))
# Resolvemos para cada valor inicial, es decir, para cada trayectoria.
t = np.linspace(0, max_time, int(250*max_time)) # Creamos un vector tiempo desde 0 hasta max_time,
# con intervalos de int(250*max_time)
x_t = np.asarray([integrate.odeint(lorenz_deriv, x0i, t)
for x0i in x0]) # resolvemos simultáneamente cada trayectoria. Empleamos un bucle for, en que resolvera
# el sistema para un primer valor, luego para el segundo y así hasta alcanzar el número
# de trayectorias introducidas. Luego almacena los resultados en una matriz.
# A lo hora de representar escogemos un color diferente para cada trayectoria.
colors = plt.cm.hot(np.linspace(0, 1, N))
# Tal y como hemos hecho en el apartado anterior, para poder representar tenemos que recuperar los datos almacenados
# en la matriz x_t. Para ello empleamos el bucle de abajo, en el que asigna a cada variable x, y, z su correspondiente
# resultado y define las características que tendrá en la gráfica, como la anchura de línea, el color...
for i in range(N):
x, y, z = x_t[i,:,:].T
lines = ax.plot(x, y, z, '-', c=colors[i])
plt.setp(lines, linewidth=2)
ax.view_init(30, angle)
plt.show()
return t, x_t
Pero vamos a observar realmente que es lo que hemos hecho. Para ello, vamos a representar las trayectorias que tanto se han nombrado:
In [7]:
t, x_t = efecto_mariposa(angle=0, N=10)
Podemos ir cambiando el valor del ángulo para girar la figura en sentido horizontal, y aumentar o disminuir el número de trayectorias modificando el valor de N en el comando anterior. Pero ya hemos visto que esta tarea la podemos implementar con un widget interactivo, por lo que... vamos a ello !!
Como habrás observado si has modificado el valor de angle y N, la representación anterior no deja muy claro el por qué se conoce a este sistema. Para observarlo más claramente tenemos que variar los parámetros del modelo. Más en concreto, para un ángulo cercano a 108 y sigma en torno a 2.5 se observa una representación que bien podrían ser consideradas las alas de una mariposa, si no lo aprecias bien, prueba a aumentar el valor de N.
In [8]:
interact(efecto_mariposa, angle=(0.,360.), N=(0,50), sigma=(0.0,50.0), rho=(0.0,50.0))
Out[8]:
Pero ¿qué sucede si variamos beta hasta alcanzar un valor cercano a 100?
Si hacemos una pequeña modificación al comando anterior, y le añadimos beta = 100, observaremos como el sistema se comporta como si estuviera dominado por un atractor tipo tórico, es decir, en forma de flotador.
In [9]:
t, x_t = efecto_mariposa(angle=0, N=10, beta = 100)
Prueba tu mismo a añadir el widget interactivo para este caso.