Processing math: 0%

domingo, 10 de julio de 2016

Integradores numéricos II: El método Leapfrog

Vamos a mostrar uno de los métodos de integración numérica más usados en sistemas conservativos. Como vimos en la entrada anterior, el integrador de Runge-Kutta no conserva la energía específica del problema de los dos cuerpos, cantidad que teóricamente se debería mantener constante en el tiempo. Se dice que este integrador no es simpléctico y en general se vuelve inapropiado para la integración del movimiento orbital a grandes periodos de tiempo. El integrador que vamos a mostrar ahora se llama Leapfrog y al ser de tipo simpléctico es capaz de preservar cantidades como la energía especifica.

Antes que nada recordemos el siguiente resultado básico del cálculo. Sea una función f(t) su expansión en serie de Taylor alrededor de una cantidad a es (a orden \mathcal{O}(t^3))

\displaystyle f(t) \approx f(a) + \left(\frac{df}{dt}\right)_{t=a}(t-a) + \frac{1}{2}\left(\frac{d^2f}{dt^2}\right)_{t=a}(t-a)^2

Ahora, consideremos una ecuación de movimiento unidimensional de la forma F(x(t),\dot{x}(t),t) = m\ddot{x}(t). Resolver esta ecuación, como sabemos, implica conocer los valores de x y \dot{x} en cada instante de tiempo. Vamos a aplicar la expansión de Taylor (alrededor del instante t) sobre la función x(t+\Delta t). De acuerdo con la ecuación anterior es claro que

x(t+\Delta t) \approx x(t) + v(t)\Delta t +\frac{1}{2}a(t)\Delta t^2

Factorizando \Delta t en los dos últimos términos de la derecha y definiendo v(t+\Delta t/2)=v(t)+\frac{1}{2}a\Delta t, la ecuación para la propagación de la posición es

x(t+\Delta t) \approx x(t) + v(t+\Delta t/2)\Delta t

Notemos que para encontrar la posición un instante de tiempo \Delta t posterior, necesitamos conocer la velocidad en un instante \Delta t/2, es decir, a mitad del intervalo de tiempo. Esa velocidad se expresa en función de la aceleración como

v(t+\Delta t/2) = v(t-\Delta t/2) + a(t)\Delta t

¿Qué pasa cuando t=0? En ese caso v(-\Delta t/2) = v(0)-a(0)\Delta t/2

Como los instantes de tiempo considerados hacen parte de un grid unidimensional de n entradas (n sería el tiempo total de integración) y ese grid es discreto, conviene reescribir las expresiones del método leapfrog en términos de un subíndice que denote el instante de tiempo considerado.

x_{n+1} = x_n + v_{n+1/2}\Delta t  \\  v_{n+1/2} = v_{n-1/2} + a_n\Delta t \\ v_{-1/2} = v_0-a_0 \Delta t/2

Las tres últimas ecuaciones constituyen el método Leapfrog, es uno de los más usados debido en gran parte a su facilidad de implementación. Lo pondremos en práctica con el mismo problema de la entrada anterior.

En primer lugar traficamos el comportamiento de la órbita en el espacio x,y y en el espacio de fases x,\dot{x}.
Figura 1.
Como podemos ver, se obtiene exactamente el mismo resultado que con el integrador odeint. Pero veamos ahora lo que sucede con la energía específica.


Figura 2.
A diferencia de lo que sucedía con el integrador odeint, en el que la energía específica aumentaba gradualmente en el tiempo, en este caso, se aprecia un comportamiento oscilante de esta misma cantidad durante el mismo intervalo de tiempo. Para propósitos prácticos, lo que interesa es conocer el valor promedio de la energía y si hacemos esto con los datos de la figura 2, encontramos que la energía en promedio se mantiene constante durante todo el tiempo de integración. Esto soporta la efectividad del Leapfrog en cuanto a la conservación de las integrales de movimiento.

#!/usr/bin/env python
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy.linalg import norm
############################################################
# Trabajamos en un sistema de unidades en el que G=1.0
G=1.0
m1=1.0
m2=0.001
mu=G*(m1+m2)
############################################################
# Vectores de condiciones iniciales y de tiempo
N=100000.0
t=np.linspace(0.0,1000.0,N)
dt=(t[len(t)-1]-t[0])/(N)
r=np.zeros((len(t),2))
v=np.zeros((len(t),2))
r[0]=[1.0, 0.0]
v[0]=[0.0, 1.2]
############################################################
# Integracion con Leap frog
for i in np.arange(1.0,len(t)):
if (i==1.0):
v[i]=v[0]-mu*r[0]/norm(r[0])**3*dt/2.0
r[i]=r[i-1]+v[i]*dt
else:
v[i]=v[i-1]-mu*r[i-1]/norm(r[i-1])**3*dt
r[i]=r[i-1]+v[i]*dt
############################################################
# Calculo de la energia especifica
Es=[]
for i,j in zip(r,v):
Es+=[0.5*norm(j)**2-mu/norm(i)]
############################################################
# Graficacion
plt.figure()
plt.plot(t,Es)
plt.xlabel('t')
plt.ylabel('Energy')
plt.title('Energy vs time')
plt.savefig('EnergyvsTime.png')
plt.figure()
plt.subplot(121)
plt.plot(r[:,0],r[:,1])
plt.xlabel('x')
plt.ylabel('y')
plt.gca().set_aspect('equal', adjustable='box') # Equal size ratio
plt.subplot(122)
plt.plot(r[:,0],v[:,0])
plt.xlabel('x')
plt.ylabel('v_x')
plt.savefig('SpaceAndPhase.png')

viernes, 24 de junio de 2016

Integradores numéricos I: El paquete SciPy

SciPy es una librería numérica de Python. Alberga una gran cantidad de rutinas de análisis numérico, desde el álgebra lineal y el cálculo hasta las ecuaciones diferenciales. Centraremos nuestra atención en la rutina odeint que nos permite resolver ecuaciones o sistemas de ecuaciones diferenciales siempre que estas sean de primer grado. Resulta que una gran cantidad de ecuaciones diferenciales que modelan fenómenos físicos son de orden dos, pero la gran mayoría de estas se pueden expresar como un conjunto de ecuaciones de primer orden a través de cambios de variable pertinentes.

Consideremos el problema gravitacional de los dos cuerpos. Cuando se describe el movimiento relativo de dos partículas sometidas a su mutua gravedad, la ecuación diferencial a resolver es

\ddot{\vec{r}}=-\frac{\mu}{r^3}\vec{r}

siendo \mu = G(m_1+m_2) y \vec{r} el vector que especifíca la posición de la partícula 2 respecto a la 1. Dado su carácter vectorial, esta ecuación es en realidad un conjunto de ecuaciones de segundo orden, una por cada grado de libertad de la partícula m_2. En principio podría pensarse que se tienen 3 grados de libertad, pero en un movimiento kepleriano, el momentum angular del sistema es una constante (en magnitud y dirección) por tanto, las dos masas se mueven sobre un mismo plano. El sistema de ecuaciones escalares es

\ddot{x}=-\frac{\mu}{r^3}x   \\  \ddot{y}=-\frac{\mu}{r^3}y

donde r=\sqrt{x^2+y^2}. Para garantizar la unicidad de la solución, debemos especificar un conjunto de condiciones iniciales, que en este caso son cuatro; las posiciones y velocidades iniciales (x_0,y_0,v_{x0},v_{y0}) de m_2 respecto a m_1. Resolver el problema, implica conocer el vector de estado Y=(x,y,v_x,v_y) en cualquier instante de tiempo.

Para solucionar numéricamente nuestro problema, el primer paso es recordar que odeint solo trabaja con ecuaciones de primer grado. Necesitamos reescribir el sistema de ecuaciones anterior apropiadamente introduciendo las variables auxiliares \dot{x}=v_x y \dot{y}=v_y, luego

\dot{x}=v_x \\ \dot{y}=v_y \\ \dot{v}_x = -\frac{\mu}{r^3}x \\ \dot{v}_y = -\frac{\mu}{r^3}y

Notemos que lo que hay a la izquierda de las igualdades anteriores, es la derivada temporal del vector de estado Y, y lo que aparece a la derecha es una función del mismo vector e implícitamente del tiempo

\frac{dY}{dt} = \text{func}(Y,t)

Una vez construida la función anterior, podemos implementarla numéricamente. Lo primero a hacer en nuestro programa es importar el módulo integrate de SciPy donde se aloja la función odeint. Esta función va a recibir necesariamente tres argumentos. El primero de ellos es la función \text{func} que como hemos mostrado, devuelve la derivada del vector de estado Y, el segundo argumento es un vector Y_0 que contiene las cuatro condiciones iniciales de nuestro problema, y el tercer argumento es un vector t, cuyas entradas son los instantes de tiempo en los cuales queremos que el integrador nos entregue la solución.

En la figura 1 hemos graficado la órbita en el espacio x,y (izquierda) y en el espacio de fases x,v_x (derecha).



Figura 1.

Vemos que con las condiciones iniciales introducidas obtenemos una órbita elíptica, con un patrón cerrado y (aparentemente) bien definido en el espacio de fases. Podría decirse entonces que el integrador ha regresado una solución correcta al problema. Pero veamos que sucede con la energía específica del sistema, que es otra cantidad conservada en el problema de los dos cuerpos

\displaystyle \epsilon = \frac{1}{2}v^2 - \frac{\mu}{r}

En la figura 2, se ve claramente que la energía específica no se conserva contrario a lo que predice la teoría. Este hecho solo puede deberse a un problema con el integrador.


Figura 2.


El asunto, es que odeint es un esquema basado en el famoso integrador Runge-Kutta. Este tipo de integradores sufre de un problema relacionado con la conservación de las propiedades de un sistema Hamiltoniano (como el problema de los dos cuerpos) en el espacio de fases. No es objeto de esta entrada discutir cuales son esas propiedades, pero lo que podemos decir es que todo integrador basado en el método de Runge-Kutta no conserva ni la energía ni las propiedades simplécticas del sistema en el espacio de fases. Esto lo hace inapropiado para cálculos de larga duración en los cuales se debe garantizar la preservación de las integrales de movimiento. En la próxima entrada estudiaremos un integrador más adecuado para tratar con este tipo de situaciones.


#!/usr/bin/env python
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy.linalg import norm
############################################################
# Trabajamos en un sistema de unidades en el que G=1.0
G=1.0
m1=1.0
m2=0.001
mu=G*(m1+m2)
############################################################
# Vectores de condiciones iniciales y de tiempo
y0=[1.0, 0.0, 0.0, 1.2]
t=np.linspace(0.0,1000.0,1000.0)
############################################################
# Definicion de la funcion que devuelve la derivada temporal
# del vector de estado 'y'
def func(y,t):
r=y[:2]
v=y[2:]
return [v[0], v[1], -mu*r[0]/norm(r)**3, -mu*r[1]/norm(r)**3]
############################################################
# Integracion con 'odeint'
solucion=odeint(func,y0,t)
r=solucion[:, :2] # vector de posicion para cada t
v=solucion[:, 2:] # vector de velocidad para cada t
############################################################
# Calculo de la energia especifica
Es=[]
for i,j in zip(r,v):
Es+=[0.5*norm(j)**2-mu/norm(i)]
############################################################
# Graficacion
plt.figure()
plt.plot(t,Es)
plt.xlabel("t")
plt.ylabel("Energy")
plt.title("Energy vs time")
plt.savefig('EnergyvsTime.png')
plt.figure()
plt.subplot(121)
plt.plot(r[:,0],r[:,1])
plt.xlabel('x')
plt.ylabel('y')
plt.gca().set_aspect('equal', adjustable='box') # Equal size ratio
plt.subplot(122)
plt.plot(r[:,0],v[:,0])
plt.xlabel('x')
plt.ylabel('v_x')
plt.savefig('SpaceAndPhase.png')


lunes, 20 de julio de 2015

Función Meshgrid y gráficos de contorno en Python

Visualizar una gráfica en 2D es relativamente sencillo, puede hacerse a punta de imaginación con algo de práctica. Pero con las funciones de tres variables la situación es más compleja. Sin embargo, siempre es posible obtener información útil a partir de una descripción bidimensional de una figura en 3D gracias a los famosos diagramas de contorno.

Imaginemos una montaña sobre la cual pasamos una serie de planos paralelos a su base a diferentes alturas. La intersección dibujará sobre cada plano la silueta que la montaña posee a esa altura. La proyección de esas siluetas sobre la base de la montaña es una diagrama de contorno. 

Formalmente, dada una función de dos variables f(x,y)=z, un diagrama de contorno se define como la curva f(x,y)=h donde h es un valor constante. En otras palabras, los diagramas de contorno son curvas donde la función adopta un valor constate. En el ejemplo de la montaña, si una persona camina a lo largo de una curva de nivel, ni asciende ni desciende. 

Aunque los diagramas de contorno son gráficos 2D, incluso en algunos casos será necesario recurrir a los dispositivos graficadores para crearlos. Vamos a dar un ejemplo. Consideremos una canica que se desliza a lo largo de una alambre parabólico (sin fricción) sometida únicamente a la acción de la gravedad. El Hamiltoniano para este sistema es

\displaystyle H(x,p) = \frac{p^2}{2(1+x^2)} + \frac{x^2}{2}

siendo p el momentum canónico asociado a la variable de posición x. Usando python y la librería matplotlib vamos a crear los diagramas de contorno de esta función. Matplotlib alberga la función contour( ), que recibe tres parámetros relacionados a través de la función H(x,p) los cuales explicaremos después. Antes de esto, notemos que nuestro gráfico de contornos estará sobre el plano x,p. Para evaluar la función H(x,p) debemos primero crear una matriz de puntos (o grid) sobre dicho plano. Para esto usamos la función meshgrid( ) de la librería numpy. Esta función recibe dos argumentos que son vectores numpy que contienen información sobre la resolución y el tamaño del grid que deseamos construir. Por ejemplo, si queremos construir una matriz de puntos que tenga un área de 20x20, los argumentos de meshgrid( ) podrían definirse así

x=np.arange(-10,11,1)
p=np.arange(-10,11,1)
view raw vectors.py hosted with ❤ by GitHub

En este caso estamos dividiendo los ejes x y p de -10 a 10 en pasos de 1 lo cual define la resolución de nuestro grid. Lo que sigue ahora es crear nuestra matriz de puntos, para lo cual hacemos

X,P=np.meshgrid(x,p)
view raw mesh.py hosted with ❤ by GitHub

Finalmente debemos evaluar la función H(x,p) en cada uno de los puntos de la matriz. Obviamente, nuestra función se evaluará en 400 puntos distintos. La función contour( ) se encarga de hacer esta evaluación automáticamente y de graficar el resultado. Los tres parámetros que recibe esta función son los tres arreglos 400-dimensionales X, Y, h  creados con ayuda de meshgrid( ). Como podemos ver, el arreglo h es precisamente la función H(x,p) evaluada en cada punto del grid. A continuación mostramos el código fuente completo y el resultado.

import numpy as np
import matplotlib.pyplot as plt
x=np.arange(-2, 2, 0.01);
p=np.arange(-2, 2, 0.01);
X,P=np.meshgrid(x,p)
h = (P**2)/(2*(1 + X**2)) + X**2/2.0
cs = plt.contour(X,P,h)
plt.clabel(cs, inline=1, fontsize=10)
plt.show()
view raw contour.py hosted with ❤ by GitHub


Notemos que hay un número que indica el valor constante de la función H(x,p) sobre cada curva. Esto se hace con la función clabel( ) de la línea diez.

martes, 4 de noviembre de 2014

El Sol y su movimiento anual

La descripción correcta del sistema solar implica abandonar la idea de que la tierra ocupa el lugar central. Aunque este hecho siempre debe estar presente en nuestra mente a veces es conveniente considerar que en efecto la Tierra es el centro. Por ejemplo, este modelo ficticio nos ayuda a responder con facilidad la siguiente pregunta: ¿por qué justo en el polo norte o en el sur el día y la noche duran 6 meses? 
Podemos considerar que durante un año, el sol se mueve a lo largo de una trayectoria (llamada eclíptica) la cual, en nuestro modelo geocéntrico, puede definirse con respecto al ecuador celeste, es decir, la proyección del ecuador de la tierra sobre una esfera imaginaria en la cual se asumen fijas todas las estrellas. Dependiendo de la época del año, el sol puede moverse arriba del ecuador celeste (en el hemisferio norte celeste) o por debajo de este (en el hemisferio sur celeste). Imaginemos ahora que estamos en el polo norte. Allí, el ecuador celeste coincide con el horizonte local del observador, por tanto, si el sol está debajo del ecuador, lo estará también del horizonte y no podremos observarlo. Esto normalmente ocurre entre septiembre y marzo. Durante los próximos seis meses, el sol estará sobre el ecuador celeste y por lo tanto sobre el horizonte, permitiendo su vista. La misma situación pero a la inversa ocurre en el polo sur.
Ahora bien, consideremos un caso más familiar: ¿cómo se mueve el Sol en el cielo de Medellín a lo largo del año? Para responder con propiedad esta pregunta, consideremos cuatro días particulares del calendario. Estos son; marzo 21 (equinoccio de primavera), Junio 21 (solsticio de verano), septiembre 22 (equinoccio de otoño) y diciembre 21 (solsticio de invierno). En los días correspondientes a los equinoccios, el Sol se encontrará justo en el ecuador celeste. 

Antes de proseguir, tengamos en cuenta que una forma de definir las coordenadas geográficas es a partir de los puntos donde la línea del ecuador celeste que surca el cielo intersecta al horizonte del observador, como lo muestra la siguiente figura.
Fig 1. Los puntos cardinales Este y Oeste pueden definirse ac
Fig 1. Los puntos cardinales Este y Oeste pueden definirse usando la línea del ecuador
celeste y su intersección con el horizonte local. 
Si el día del equinoccio, el Sol se mueve a lo largo del ecuador celeste y este define las coordenadas oriente y occidente, entonces es claro que solo en esos días, el sol saldrá exactamente por el Este y se pondrá justo por el Oeste.
Ahora estudiemos los solsticios. Como ya se dijo, existen días en los cuales la posición del sol coincide con el ecuador celeste y otros en los que estará sobre o por debajo de él. Este último caso corresponde a los solsticios. Para un observador en el hemisferio norte, en el día del solsticio de verano, el 21 de junio aproximadamente, el Sol estará sobre el ecuador celeste a exactamente 23,5º. Por esta razón se dice, que en el día del solsticio de verano, el Sol alcanza su punto más alto en el cielo. Lo opuesto sucede el 21 de diciembre, en el solsticio de invierno, donde el Sol alcanza una altura mínima.


Fig 2. Posición del Sol el día del solsticio de invierno
en el hemisferio norte.

Ahora, hagamos un ejercicio mental. Imaginemos una línea paralela al ecuador que pase justo sobre el cenit de un observador situado a una latitud \delta. Si \delta>0º , es decir, si estamos en el hemisferio norte, con respecto a esta línea imaginaria el Sol podrá moverse 23,5º + \delta hacia el sur y 23,5º - \delta hacia el norte. En este caso, la condición necesaria para que el Sol pase por el cenit es 
23,5º - |\delta|\geq 0º 
Si estamos en el hemisferio sur, con respecto a una paralela al ecuador que pase por el cenit, el Sol podrá moverse 23,5º - \delta hacia el sur y 23,5º + \delta hacia el norte. En este caso, la condición de paso por el cenit es exactamente la misma. El valor absoluto de la latitud deriva del hecho de que para el hemisferio sur \delta<0 y para el norte \delta>0.

Contrario a lo que comúnmente se cree, el Sol no se encuentra justo sobre nuestras cabezas al medio día durante todo el año. Solo en días particulares ocurrirá este fenómeno que es conocido como el día sin sombra. Este evento solo es observable para latitudes entre los 23,5º N y 23,5º S. En la ciudad de Medellín, pudimos presenciarlo el pasado 6 de septiembre. 



Fig 3. Cortesía: Planetario de Medellín 



lunes, 24 de febrero de 2014

Testing MathJax II

Vamos a implementar MathJax a nuestras entradas. En primer lugar, luego de loguearnos en blogger y seleccionar el blog con el que deseamos trabajar, damos clic en Plantilla, a la izquierda de la pantalla. Luego, clic en Editar HTML lo que desplegará el código fuente del sitio. Vamos a modificar un poco este código. Antes que nada es importante advertir sobre el cuidado que debe tener al momento de trabajar con el código fuente. Borrar o modificar un elemento por insignificante que parezca puede resultar fatal si no sabe lo que está haciendo. Una vez dentro del código buscamos la etiqueta <head> la cuál normalmente se encuentra en los primeros renglones. Justo debajo de ella, pegamos en siguiente código y guardamos los cambios.
Si su plantilla es simple, lo que ha hecho hasta el momento bastará. Si es dinámica (como la de este blog) debe agregar el siguiente código (usando el editor HTML) al final de cada post que utilice \LaTeX.
<script src='http://cdn.mathjax.org/mathjax/latest/MathJax.js' type='text/javascript'>
MathJax.Hub.Config({
HTML: [&quot;input/TeX&quot;,&quot;output/HTML-CSS&quot;],
TeX: { extensions: [&quot;AMSmath.js&quot;,&quot;AMSsymbols.js&quot;],
equationNumbers: { autoNumber: &quot;AMS&quot; } },
extensions: [&quot;tex2jax.js&quot;],
jax: [&quot;input/TeX&quot;,&quot;output/HTML-CSS&quot;],
tex2jax: { inlineMath: [ [&#39;&#39;,&#39;&#39;], [&quot;\\(&quot;,&quot;\\)&quot;] ],
displayMath: [ [&#39;&#39;,&#39;&#39;], [&quot;\\[&quot;,&quot;\\]&quot;] ],
processEscapes: true },
&quot;HTML-CSS&quot;: { availableFonts: [&quot;TeX&quot;],
linebreaks: { automatic: true } }
});
</script>
view raw mathjax.html hosted with ❤ by GitHub
Eso es todo, ahora estamos en capacidad de escribir expresiones como \log_{10} \left[\frac{I(r)}{I_e}\right]=-3,3307\left[\left(\frac{r}{r_e}\right)^{1/4}-1\right] Créditos a Paulo Loma por su investigación sobre este tema.

sábado, 22 de febrero de 2014

Testing MathJax


Como mencioné en la entrada anterior, espero que este blog me ayude a compartir algunas de las cosas que me apasionan. Temas relacionados con las ciencias físicas son de particular interés para mí. Si eventualmente debo escribir alguna expresión matemática, me gustaría que tuviera una buena calidad en cuanto a tipografía. Ni Blogger, ni Wordpress ni ningún otro servicio de "blogging" gratuito ofrece por si solo la posibilidad de anidar tipografía estilo \LaTeX en las entradas. En todos los casos es necesario incluir bibliotecas especiales que permitan hacer esto. MathJax es un excelente ejemplo. Se trata de una librería basada en JavaScript que integra notación matemática en los navegadores web por medio de lenguajes de marcado como \LaTeX, MathML y ASCIIMathML. Es decir, al escribir en nuestro panel de edición la expresión $ \sin²(\theta) + \cos²(\theta) = 1 $, el usuario verá lo siguiente:

\sin^2(\theta) + \cos^2(\theta) = 1

Lograr lo anterior requirió un par de horas de "Googling", así que la siguiente entrada estará dedicada a mostrar paso a paso el procedimiento.

viernes, 21 de febrero de 2014

Hello world!


Desde hace tiempo quería tener un espacio que me permitiera compartir algunas de las cosas que considero interesantes y también otras que no lo son tanto pero que aún así debo hacerlas. Esto de los blogs se ha vuelto muy de moda últimamente y no quería resagarme respecto a las nuevas corrientes, más aún si considero la utilidad de esta herramienta al momento de llevar un registro de algunas de mis tareas cotidianas que muchas veces, sea por pereza o por falta de tiempo, dejo inconclusas. En realidad, creo que escribo este blog pensando más en mí que en los demás, así que al menos por ahora, no tendrá mayor importancia si lo que publico aquí llega o no a más personas.

Si luego de leer lo anterior aún desea seguir mis actualizaciones o revisar este sitio con frecuencia, ¡sea usted bienvenido!