Un poco de mecánica celeste

Los cálculos, como es ya habitual, se los dejaremos a SymPy, pero esta vez usaremos IPython como interfaz de SymPy. IPython permite trabajar con notebooks, que son documentos que se visualizan con un navegador y permiten mezclar código Python (u otros lenguajes), con código html y LateX. La presente notebook se aloja en el siguiente enlace.

Los mismos cálculos se efectuaron con SAGE, y se pueden encontrar en el script 2cuerpos.sage dentro de la carpeta scripts en el repositorio de GitHub que mantiene materiales de este curso. Ver el siguiente enlace

Supondremos que el movimiento del planeta se realiza en un plano. Esta suposición no entraña una pérdida de generalidad, dado que es posible demostrar que es cierta. Esta demostración es sencilla, pero nos alejaría de nuestro objetivo, por tanto la omitiremos. Vamos a usar coordenadas polares $(r,\theta)$ y los versores $\overrightarrow{u}_r:=(\cos\theta,\sin\theta)$ y $\overrightarrow{u}_{\theta}:=(-\sin\theta, \cos\theta)$. Notar que $\overrightarrow{u}_r \perp \overrightarrow{u}_{\theta}$ y por consiguiente $\mathcal{B}:=\{\overrightarrow{u}_r , \overrightarrow{u}_{\theta}\}$ forma una base del espacio euclideano 2-dimensional. Usaremos este hecho para representar distintos vectores como combinación lineal de vectores de la base.



In [1]:

    
from sympy import *
t,mu=symbols('t,mu')
x=Function('x')(t)
y=Function('y')(t)
r=Function('r')(t)
theta=Function('theta')(t)
u_r=Matrix([cos(theta),sin(theta)])
init_printing()
u_r









    Out[1]:





$$\left[\begin{matrix}2 \cos{\left (\theta{\left (t \right )} \right )}\\2 \sin{\left (\theta{\left (t \right )} \right )}\end{matrix}\right]$$



In [4]:

    
u_theta=Matrix([-sin(theta),cos(theta)])
u_theta









    Out[4]:





$$\left[\begin{matrix}- \sin{\left (\theta{\left (t \right )} \right )}\\\cos{\left (\theta{\left (t \right )} \right )}\end{matrix}\right]$$

Vamos a representar el vector aceleración $\overrightarrow{a}=\frac{d^2\overrightarrow{u}}{dt^2}$



In [6]:

    
pos=r*u_r
a1,a2=symbols('a1,a2')
sol=solve(a1*u_r+a2*u_theta-pos.diff(t,2),[a1,a2])
sol









    Out[6]:





$$\begin{Bmatrix}a_{1} : - r{\left (t \right )} \left(\frac{d}{d t} \theta{\left (t \right )}\right)^{2} + \frac{d^{2}}{d t^{2}}  r{\left (t \right )}, & a_{2} : r{\left (t \right )} \frac{d^{2}}{d t^{2}}  \theta{\left (t \right )} + 2 \frac{d}{d t} r{\left (t \right )} \frac{d}{d t} \theta{\left (t \right )}\end{Bmatrix}$$



In [7]:

    
sol[a1]









    Out[7]:





$$- r{\left (t \right )} \left(\frac{d}{d t} \theta{\left (t \right )}\right)^{2} + \frac{d^{2}}{d t^{2}}  r{\left (t \right )}$$



In [8]:

    
sol[a2]









    Out[8]:





$$r{\left (t \right )} \frac{d^{2}}{d t^{2}}  \theta{\left (t \right )} + 2 \frac{d}{d t} r{\left (t \right )} \frac{d}{d t} \theta{\left (t \right )}$$

El vector aceleración debe ser igual a la fuerza por unidad de masa $-\mu \overrightarrow{r}/r^3$. Notemos que esta fuerza es central, es decir tiene componente nula respecto al vector $\overrightarrow{u}_{\theta}$. Por consiguiente se debe satisfacer que

$$\frac{1}{r}\frac{d}{dt} \left( r^2\dot{\theta} \right)=0\Longleftrightarrow \exists h\in\mathbb{R}: \boxed{ r^2\dot{\theta}=h}. $$

Hemos derivado la Segunda Ley de Kepler: El radio vector barre áreas iguales en tiempos iguales.

En la dirección radial $\overrightarrow{u}_r$ la componente de la fuerza es $-\mu/r^2$. Es decir se satisface la ecuación $$ \ddot{r}-r\dot{\theta}^2=-\frac{\mu}{r^2} $$

Notar que esta ecuación entraña dos incognitas $r$ y $\theta$, pero $\dot{\theta}$ puede ser remplazado por $h/r^2$ por la segunda Ley de Kepler. Declaremos la variable $h$ que juega un rol importante y reemplacemos $\dot{\theta}$ en la ecuación



In [9]:

    
h=symbols('h')
ed=(sol[a1]).subs(theta.diff(t),h/r**2)
ed









    Out[9]:





$$- \frac{h^{2}}{r^{3}{\left (t \right )}} + \frac{d^{2}}{d t^{2}}  r{\left (t \right )}$$

Conseguimos una ecuación no lineal de segundo orden para $r$. De los métodos que hemos visto, ninguno se aplica a esta ecuación. El truco mágico consiste en considerar la nueva variable dependiente $z=1/r$ y la nueva variable independiente $\theta$.



In [10]:

    
z=Function('z')(theta)
r=1/z
ed2=r.diff(t,2)+mu/r**2-h**2/r**3
ed2









    Out[10]:





$$- h^{2} z^{3}{\left (\theta{\left (t \right )} \right )} + \mu z^{2}{\left (\theta{\left (t \right )} \right )} + \frac{1}{z^{2}{\left (\theta{\left (t \right )} \right )}} \left(- \left(\frac{d}{d t} \theta{\left (t \right )}\right)^{2} \left. \frac{d^{2}}{d \xi_{1}^{2}}  z{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=\theta{\left (t \right )} }} - \frac{d^{2}}{d t^{2}}  \theta{\left (t \right )} \left. \frac{d}{d \xi_{1}} z{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=\theta{\left (t \right )} }} + \frac{2 \left(\frac{d}{d t} \theta{\left (t \right )}\right)^{2}}{z{\left (\theta{\left (t \right )} \right )}} \left(\left. \frac{d}{d \xi_{1}} z{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=\theta{\left (t \right )} }}\right)^{2}\right)$$

En la ecuación resultante, nuevamente aparece $\dot{\theta}$ y además ahora aparece $\ddot{\theta}$. Tenemos que reemplazar $\dot{\theta}$ por $hz^2$ y $\ddot{\theta}$ por $\tfrac{d}{dt}hz^2$.



In [11]:

    
theta2diff=(h*z**2).diff(t).subs(theta.diff(t),h*z**2)
ed3=ed2.subs(theta.diff(t,2),theta2diff)
ed3=ed3.subs(theta.diff(t),h*z**2)
ed3









    Out[11]:





$$- h^{2} z^{3}{\left (\theta{\left (t \right )} \right )} - h^{2} z^{2}{\left (\theta{\left (t \right )} \right )} \left. \frac{d^{2}}{d \xi_{1}^{2}}  z{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=\theta{\left (t \right )} }} + \mu z^{2}{\left (\theta{\left (t \right )} \right )}$$



In [12]:

    
(ed3/h**2/z**2).expand()









    Out[12]:





$$- z{\left (\theta{\left (t \right )} \right )} - \left. \frac{d^{2}}{d \xi_{1}^{2}}  z{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=\theta{\left (t \right )} }} + \frac{\mu}{h^{2}}$$



In [ ]: