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')

view raw 2body-NoSymplectic.py hosted with ❤ by GitHub