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. |
![]() |
Figura 2. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/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') |