Un método para graficar órbitas y trayectorias en la teoría Newtoniana de la Gravitación

[Richard R Richard]

Voy a trabajar con un caso sencillo, el problema de dos cuerpos en 3 dimensiones, donde el objeto más masivo tiene centro en el origen de coordenadas, pero es extensible a $n$ cuerpos a modificaciones en el origen de coordenadas, incluso poniéndolas en el CM del sistema.

Un diagrama de flujo típico de programación en ordenadores puede ser el siguiente

Donde empezamos definiendo constantes como $G$ que es la constante de gravitación universal y $M$ se interpreta como la masa del objeto celeste más masivo,( planeta, estrella etc) que crea el gravedad o algunas otra útiles como Pi si no vienen predefinidas en el el lenguaje.

Para definir una trayectoria para una partícula de prueba hay que dotarla de los parámetros iniciales de su movimiento. es decir hay que definir su posición y su velocidad para un instante inicial dado en coordenadas esféricas o bien, más fácil aún, en cartesianas . Adoptare coordenadas cartesianas.por lo que podemos convertir desde sistemas de coordenadas esféricas

$$\overline x(t_0)=(r_0,\theta_0,\phi_0)$$

a $\overline x=(x_0,y_0,z_0)$

usando una conversión de coordenadas

$$\left[\begin{aligned}<br /> x & = r\sin\theta\cos\phi \\<br /> y & = r\sin\theta\sin\phi \\<br /> z & = r\cos\theta \end{aligned}\right ]$$

y recalcular velocidades

$$\overline v=(v_{x0},v_{y0},v_{z0})=( \dot x_{(t_0)},\dot y_{(t_0)},\dot y_{(t_0)})=(u_1,u_2,u_3)$$

Donde $\dfrac{\partial x_i}{\partial t }=u_i =\dot x_i$ y $\dfrac{\partial^2 x_i}{\partial t^2 }=a_i =\ddot x_i$


Definimos la posición del objeto masivo a tiempo $\vec x_M(t)=(0,0,0)\quad \forall t>0$ y establecemos que no variara durante el tiempo que simulemos el problema.

El número de iteraciones no puede ser tan grande como queremos, ya que los ordenadores introducen un error de truncamiento en la última cifra de los resultado presentados, luego el número de iteraciones multiplicado por el valor del error de truncamiento , debe permanecer por debajo del error aceptable en el resultado de la medición.

Lo que sigue es el cálculo de la aceleración de la partícula de prueba en el campo gravitatorio

$$\dfrac{\partial ^2 r}{\partial t^2}=-\dfrac{GM}{r^2}$$

Luego podemos hacer proyección en cualquiera de los ejes usando

$\cos \theta=\dfrac{x}{|r|}$ , $\sin \theta=\dfrac{y}{|r|}$ , $\sin\phi=\dfrac{z}{|r|}$

Por lo que la proyección de la aceleración en cada eje cartesiano puede ser calculada como

$$\dfrac{\partial ^2 x}{\partial t^2}=-\dfrac{GM}{r^2}\cos \theta=-\dfrac{GMx}{|r|^3}=-\dfrac{GMx}{(x^2+y^2+z^2)^{\frac32}}$$

$$\dfrac{\partial ^2 y}{\partial t^2}=-\dfrac{GM}{r^2}\sin \theta=-\dfrac{GMy}{|r|^3}=-\dfrac{GMy}{(x^2+y^2+z^2)^{\frac32}}$$

$$\dfrac{\partial ^2 z}{\partial t^2}=-\dfrac{GM}{r^2}\sin \phi=-\dfrac{GMz}{|r|^3}=-\dfrac{GMz}{(x^2+y^2+z^2)^{\frac32}}$$

La velocidad de la partícula en cada iteración será evaluada incrementando el tiempo en $\Delta t =\dfrac{T_{total}}{paso}$

resultando

$$\left[\begin{aligned} t(i) &= i \cdot \Delta t \\<br /> u_x(t_i) &= u(t_{i-1})_x+\dfrac{a_x(t_i)-a_x(t_{i-1})}{2}\Delta t \\<br /> u_y(t_i) &= u(t_{i-1})_y+\dfrac{a_y(t_i)-a_y(t_{i-1})}{2}\Delta t \\<br /> u_z(t_i) &= u(t_{i-1})_z+\dfrac{a_z(t_i)-a_z(t_{i-1})}{2}\Delta t \end{aligned}\right]$$

como verán uso la aceleracion media y no la aceleración a un tiempo anterior pues he visto mejores resultados así

Y para la posición

$$\left[\begin{aligned}x(t_i) &= x(t_{i-1})+u_x(t_{i-1})\Delta t+\dfrac{a_x(t_i)}{2}\Delta^2 t \\<br /> y(t_i) &= y(t_{i-1})+u_y(t_{i-1})\Delta t+\dfrac{a_y(t_i)}{2}\Delta^2 t \\<br /> z(t_i) &= z(t_{i-1})+u_z(t_{i-1})\Delta t+\dfrac{a_z(t_i)}{2}\Delta^2 t \end{aligned}\right]$$

Luego de registrar esos valores en la tabla de posiciones, en función del tiempo $t_i$ pasamos a actualizar variables con $i=i+1$ y preguntándonos si hemos llegado al tiempo total $T_{total}$.

Luego queda reproducir en escala de tiempo real. T=12 o 24 fps o acelerada donde el salto entre frame y frame es proporcional a otro salto de tiempo. Esto ira dejando dando idea de la trayectoria de punto a punto, pero se puede dejar unidos por rectas, (o curvas usando métodos mas complejos) para dejar la órbita descrita de forma mas aproximada.

La representación 2D en pantalla proveniente de datos en 3D , se aborda convirtiendo, coordenadas a pixeles mediante una matriz de conversión dependiendo del tipo de proyección que busquemos,


Proyección central con un punto de fuga.
Proyección paralela

• isométrica.
• Dimétrica
• Trimétrica

Proyección oblicua: Persp. caballera.
Proyección oblicua:Perspectiva militar.

Con las ultimas dos no se obtienen resultados muy realistas al estilo de los videojuegos.

[arivasm]

Hace años trabajé en este tema. Recuerdo que el mejor procedimiento para la integración era recurrir a los métodos de Runge-Kutta (quizá ahora ya haya cosas incluso mejores). Puedes chequear la calidad del método iterativo que emplees revisando el cumplimiento de las leyes de conservación. Si no me equivoco, el algoritmo que propones tiene el problema de que dichas magnitudes (la energía, por ejemplo) se van apartando cada vez más rápido del valor original, mientras que con un Runge-Kutta de orden mínimamente elevado (creo que yo usaba el 4, pero también me suena haber implementado órdenes mayores con buenos resultados) oscila, aunque la amplitud de las oscilaciones va aumentando con las iteraciones.

[Richard R Richard]

Si claro arivasm para mejorar algunos cálculos que hice en su momento, usé Runge- Kutta, se puede ver en Análisis Numérico de Richard L Burden y J Douglas Faires pag 243 ISBN 968-7270-09-08, ese libro me ha salvado un dolor de cabeza muchísimas veces.