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

view raw 2Body-Symplectic.py hosted with ❤ by GitHub