simulación de N cuerpos

En los métodos de árbol , como una simulación de Barnes-Hut , se suele usar un octree para dividir el volumen en celdas cúbicas y solo las interacciones entre partículas de celdas cercanas necesitan ser tratadas individualmente; las partículas en celdas distantes pueden tratarse colectivamente como una sola partícula grande centrada en el centro de masa de la celda distante (o como una expansión multipolar de bajo orden ). Esto puede reducir drásticamente la cantidad de interacciones de pares de partículas que deben calcularse. Para evitar que la simulación se sature al calcular interacciones partícula-partícula, las celdas deben refinarse a celdas más pequeñas en las partes más densas de la simulación que contienen muchas partículas por celda. Para simulaciones donde las partículas no están distribuidas uniformemente, los métodos de descomposición de pares bien separados de Callahan y Kosaraju producen un tiempo óptimo O( n  log n ) por iteración con dimensión fija. 

Otra posibilidad es el método de malla de partículas, en el que el espacio se discretiza en una malla y, para calcular el potencial gravitatorio , se supone que las partículas se dividen entre los vértices circundantes de 2x2 de la malla. La energía potencial Φ se puede obtener mediante la ecuación de Poisson.

$${\displaystyle \nabla ^{2}\Phi =4\pi G{\rho },\,}$$

donde G es la constante de Newton y${\displaystyle \rho }$es la densidad (número de partículas en los puntos de la malla). La transformada rápida de Fourier puede resolver esto de manera eficiente pasando al dominio de la frecuencia donde la ecuación de Poisson tiene la forma simple

$${\displaystyle {\sombrero {\Phi }}=-4\pi G{\frac {\sombrero {\rho }}{k^{2}}},}$$

dónde${\displaystyle {\vec {k}}}$es el número de onda comóvil y los sombreros denotan transformadas de Fourier. Dado que${\displaystyle {\vec {g}}=-{\vec {\nabla }}\Phi }$, ahora se puede encontrar el campo gravitatorio multiplicando por${\displaystyle -i{\vec {k}}}$y calculando la transformada inversa de Fourier (o calculando la transformada inversa y luego utilizando algún otro método). Dado que este método está limitado por el tamaño de la malla, en la práctica se utiliza una malla más pequeña u otra técnica (como la combinación con un árbol o un algoritmo simple de partículas-partículas) para calcular las fuerzas a pequeña escala. A veces se utiliza una malla adaptativa, en la que las celdas de la malla son mucho más pequeñas en las regiones más densas de la simulación.

Se utilizan varios algoritmos diferentes de perturbación gravitacional para obtener estimaciones bastante precisas de la trayectoria de los objetos en el Sistema Solar .

A menudo se opta por colocar un satélite en una órbita fija . La trayectoria de un satélite que orbita cerca de la Tierra se puede modelar con precisión partiendo de una órbita elíptica de dos cuerpos alrededor del centro de la Tierra, añadiendo pequeñas correcciones debidas al achatamiento terrestre , la atracción gravitatoria del Sol y la Luna, la resistencia atmosférica, etc. Es posible encontrar una órbita fija sin calcular la trayectoria real del satélite.

La trayectoria de un planeta pequeño, un cometa o una nave espacial de largo alcance a menudo se puede modelar con precisión partiendo de la órbita elíptica de dos cuerpos alrededor del Sol y añadiendo pequeñas correcciones debidas a la atracción gravitatoria de los planetas más grandes en sus órbitas conocidas.

Algunas características de las trayectorias a largo plazo de un sistema de partículas pueden calcularse directamente. No es necesario calcular la trayectoria real de ninguna partícula en particular como paso intermedio. Dichas características incluyen la estabilidad de Lyapunov , el tiempo de Lyapunov , diversas mediciones de la teoría ergódica , etc.

El suavizado es un truco numérico utilizado en las técnicas de N cuerpos para evitar divergencias numéricas cuando una partícula se acerca demasiado a otra (y la fuerza tiende al infinito). Esto se obtiene modificando el potencial gravitatorio regularizado de cada partícula como

$${\displaystyle \Phi =-{\frac {1}{\sqrt {r^{2}+\epsilon ^{2}}}},}$$

(en lugar de 1/r) donde${\displaystyle \epsilon }$es el parámetro de suavizado. El valor del parámetro de suavizado debe establecerse lo suficientemente pequeño para que las simulaciones sean realistas.

La implementación más simple de simulaciones de N cuerpos donde${\textstyle n\geq 3}$es una propagación ingenua de cuerpos en órbita; ingenua implica que las únicas fuerzas que actúan sobre los cuerpos en órbita son las fuerzas gravitatorias que ejercen entre sí. En lenguajes de programación orientados a objetos , como C++ , cierto código repetitivo es útil para establecer las estructuras matemáticas fundamentales, así como los contenedores de datos necesarios para la propagación; a saber, vectores de estado y, por lo tanto, vectores , y algún objeto fundamental que contenga estos datos, así como la masa de un cuerpo en órbita. Este método también es aplicable a otros tipos de simulaciones de N cuerpos; una simulación de masas puntuales con cargas utilizaría un método similar, sin embargo, la fuerza se debería a la atracción o repulsión por interacción de campos eléctricos. En cualquier caso, la aceleración de una partícula es el resultado de la suma de los vectores de fuerza, dividida por la masa de la partícula:

$${\displaystyle {\vec {a}}={\frac {1}{m}}\sum {\vec {F}}}$$

Un ejemplo de un método programáticamente estable y escalable para contener datos cinemáticos de una partícula es el uso de matrices de longitud fija, lo que en código optimizado permite una fácil asignación de memoria y predicción de los recursos consumidos; como se puede ver en el siguiente código C++:

struct Vector3{doble e [ 3 ] = { 0 };Vector3 () {}~ Vector3 () {}Vector3 en línea ( doble e0 , doble e1 , doble e2 ){esto -> e [ 0 ] = e0 ;esto -> e [ 1 ] = e1 ;esto -> e [ 2 ] = e2 ;}};struct OrbitalEntity{doble e [ 7 ] = { 0 };OrbitalEntity () {}~ OrbitalEntity () {}inline OrbitalEntity ( double e0 , double e1 , double e2 , double e3 , double e4 , double e5 , double e6 ){esto -> e [ 0 ] = e0 ;esto -> e [ 1 ] = e1 ;esto -> e [ 2 ] = e2 ;esto -> e [ 3 ] = e3 ;esto -> e [ 4 ] = e4 ;esto -> e [ 5 ] = e5 ;esto -> e [ 6 ] = e6 ;}};

Tenga en cuenta que OrbitalEntitycontiene suficiente espacio para un vector de estado, donde:

• ${\textstyle e_{0}=x}$, la proyección del vector de posición de los objetos en el espacio cartesiano a lo largo de${\displaystyle \left[1\;0\;0\right]}$
• ${\textstyle e_{1}=y}$, la proyección del vector de posición de los objetos en el espacio cartesiano a lo largo de${\displaystyle \left[0\;1\;0\right]}$
• ${\textstyle e_{2}=z}$, la proyección del vector de posición de los objetos en el espacio cartesiano a lo largo de${\displaystyle \left[0\;0\;1\right]}$
• ${\textstyle e_{3}={\dot {x}}}$, la proyección del vector de velocidad del objeto en el espacio cartesiano a lo largo de${\displaystyle \left[1\;0\;0\right]}$
• ${\textstyle e_{4}={\dot {y}}}$, la proyección del vector de velocidad del objeto en el espacio cartesiano a lo largo de${\displaystyle \left[0\;1\;0\right]}$
• ${\textstyle e_{5}={\dot {z}}}$, la proyección del vector de velocidad del objeto en el espacio cartesiano a lo largo de${\displaystyle \left[0\;0\;1\right]}$

Además, OrbitalEntitycontiene espacio suficiente para un valor de masa.

Comúnmente, las simulaciones de N cuerpos serán sistemas basados ​​en algún tipo de ecuaciones de movimiento ; de estas, la mayoría dependerán de alguna configuración inicial para "sembrar" la simulación. En sistemas como aquellos que dependen de algún potencial gravitatorio o eléctrico, la fuerza sobre una entidad de simulación es independiente de su velocidad. Por lo tanto, para sembrar las fuerzas de la simulación, solo se necesitan posiciones iniciales, pero esto no permitirá la propagación: se requieren velocidades iniciales. Consideremos un planeta orbitando una estrella: no tiene movimiento, pero está sujeto a la atracción gravitatoria de su estrella anfitriona. A medida que avanza el tiempo y se agregan pasos de tiempo , ganará velocidad de acuerdo con su aceleración. Para un instante dado en el tiempo,${\displaystyle t_{n}}$, la aceleración resultante de un cuerpo debido a sus masas vecinas es independiente de su velocidad; sin embargo, para el paso de tiempo${\displaystyle t_{n+1}}$, el cambio de posición resultante es significativamente diferente debido a la dependencia inherente de la propagación con respecto a la velocidad. En los mecanismos de propagación básicos, como el método de Euler simpléctico que se utilizará más adelante, la posición de un objeto en${\displaystyle t_{n+1}}$depende únicamente de su velocidad en${\displaystyle t_{n}}$, ya que el cambio de posición se calcula mediante

$${\displaystyle {\vec {r}}_{t_{n+1}}={\vec {r}}_{t_{n}}+{\vec {v}}_{t_{n}}\cdot \Delta t}$$

Sin aceleración,${\textstyle {\vec {v}}_{t_{n}}}$es estático, sin embargo, desde la perspectiva de un observador que solo ve la posición, se necesitarán dos pasos de tiempo para ver un cambio en la velocidad.

Se puede realizar una simulación similar a la del sistema solar calculando las distancias promedio de masas puntuales equivalentes a planetas desde una estrella central. Para simplificar el código, se utilizará un enfoque no riguroso basado en semiejes mayores y velocidades medias. Es necesario reservar espacio de memoria para estos cuerpos antes de configurarlos; para permitir la escalabilidad, se puede usar el comando malloc .

OrbitalEntity * orbital_entities = malloc ( sizeof ( OrbitalEntity ) * ( 9 + N_ASTEROIDS ));orbital_entities [ 0 ] = { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 1.989e30 }; // una estrella similar al solorbital_entities [ 1 ] = { 57.909e9 , 0.0 , 0.0 , 0.0 , 47.36e3 , 0.0 , 0.33011e24 }; // un planeta similar a Mercurioorbital_entities [ 2 ] = { 108.209e9 , 0.0 , 0.0 , 0.0 , 35.02e3 , 0.0 , 4.8675e24 }; // un planeta similar a Venusorbital_entities [ 3 ] = { 149.596e9 , 0.0 , 0.0 , 0.0 , 29.78e3 , 0.0 , 5.9724e24 }; // un planeta similar a la Tierraorbital_entities [ 4 ] = { 227.923e9 , 0.0 , 0.0 , 0.0 , 24.07e3 , 0.0 , 0.64171e24 }; // un planeta similar a Marteorbital_entities [ 5 ] = { 778.570e9 , 0.0 , 0.0 , 0.0 , 13e3 , 0.0 , 1898.19e24 }; // un planeta similar a Júpiterorbital_entities [ 6 ] = { 1433.529e9 , 0.0 , 0.0 , 0.0 , 9.68e3 , 0.0 , 568.34e24 }; // un planeta similar a Saturnoorbital_entities [ 7 ] = { 2872.463e9 , 0.0 , 0.0 , 0.0 , 6.80e3 , 0.0 , 86.813e24 }; // un planeta similar a Uranoorbital_entities [ 8 ] = { 4495.060e9 , 0.0 , 0.0 , 0.0 , 5.43e3 , 0.0 , 102.413e24 }; // un planeta similar a Neptuno

donde N_ASTEROIDSes una variable que permanecerá en 0 temporalmente, pero permite la inclusión futura de un número significativo de asteroides, a discreción del usuario. Un paso crítico para la configuración de las simulaciones es establecer los rangos de tiempo de la simulación,${\displaystyle t_{0}}$a${\displaystyle t_{\text{fin}}}$así como el paso de tiempo incremental${\displaystyle dt}$lo que hará avanzar la simulación:

doble t_0 = 0 ;doble t = t_0 ;doble dt = 86400 ;double t_end = 86400 * 365 * 10 ; // aproximadamente una década en segundosdoble BIG_G = 6,67e-11 ; // constante gravitacional

Las posiciones y velocidades establecidas anteriormente se interpretan como correctas para${\displaystyle t=t_{0}}$.

El alcance de una simulación lógicamente sería para el período donde${\displaystyle t_{0}\leq t<t_{\text{fin}}}$.

Una simulación completa puede constar de innumerables pasos de tiempo. En el nivel elemental, cada paso de tiempo implica calcular, para cada cuerpo:

• las fuerzas que actúan sobre el cuerpo;
• aceleración (${\displaystyle {\vec {a}}}$);
• velocidad (${\displaystyle {\vec {v}}_{n}={\vec {v}}_{n-1}+{\vec {a}}_{n}\cdot \Delta t}$);
• y posición (${\displaystyle {\vec {r}}_{n+1}={\vec {r}}_{n}+{\vec {v}}_{n}\cdot \Delta t}$).

Lo anterior se puede implementar de forma bastante sencilla con un bucle while que continúa mientras${\displaystyle t}$existe en el rango mencionado:

mientras ( t < t_end ) { para ( size_t m1_idx = 0 ; m1_idx < 9 + N_ASTEROIDS ; m1_idx ++ ) { Vector3 a_g = { 0 , 0 , 0 };para ( tamaño_t m2_idx = 0 ; m2_idx < 9 + N_ASTEROIDS ; m2_idx ++ ) { si ( m2_idx != m1_idx ) { Vector3 r_vector ;r_vector . e [ 0 ] = orbital_entities [ m1_idx ]. e [ 0 ] - orbital_entities [ m2_idx ]. e [ 0 ]; r_vector . e [ 1 ] = orbital_entities [ m1_idx ]. e [ 1 ] - orbital_entities [ m2_idx ]. e [ 1 ]; r_vector . e [ 2 ] = orbital_entities [ m1_idx ]. e [ 2 ] - orbital_entities [ m2_idx ]. e [ 2 ];double r_mag = sqrt ( r_vector . e [ 0 ] * r_vector . e [ 0 ] + r_vector . e [ 1 ] * r_vector . e [ 1 ] + r_vector . e [ 2 ] * r_vector . e [ 2 ]);doble aceleración = -1.0 * BIG_G * ( orbital_entities [ m2_idx ]. e [ 6 ]) / pow ( r_mag , 2.0 );Vector3 r_unit_vector = { r_vector . e [ 0 ] / r_mag , r_vector . e [ 1 ] / r_mag , r_vector . e [ 2 ] / r_mag };a_g . e [ 0 ] += aceleración * r_unit_vector . e [ 0 ]; a_g . e [ 1 ] += aceleración * r_unit_vector . e [ 1 ]; a_g . e [ 2 ] += aceleración * r_unit_vector . e [ 2 ]; } }orbital_entities [ m1_idx ] .e [ 3 ] + = a_g.e [ 0 ] * dt ; orbital_entities [ m1_idx ] .e [ 4 ] + = a_g.e [ 1 ] * dt ; orbital_entities [ m1_idx ] .e [ 5 ] + = a_g.e [ 2 ] * dt ; }​​​for ( size_t entity_idx = 0 ; entity_idx < 9 + N_ASTEROIDS ; entity_idx ++ ) { orbital_entities [ entity_idx ]. e [ 0 ] += orbital_entities [ entity_idx ]. e [ 3 ] * dt ; orbital_entities [ entity_idx ]. e [ 1 ] += orbital_entities [ entity_idx ]. e [ 4 ] * dt ; orbital_entities [ entity_idx ]. e [ 2 ] += orbital_entities [ entity_idx ]. e [ 5 ] * dt ; } t += dt ; }

Centrándonos en los cuatro planetas rocosos interiores de la simulación, las trayectorias resultantes de la propagación anterior se muestran a continuación:

• von Hoerner, Sebastián (1960). "Die numerische Integration des n-Körper-Problemes für Sternhaufen. I". Zeitschrift für Astrophysik (en alemán). 50 : 184. Código Bib : 1960ZA.....50..184V .
• von Hoerner, Sebastián (1963). "Die numerische Integration des n -Körper-Problemes für Sternhaufen. II". Zeitschrift für Astrophysik (en alemán). 57 : 47. Código Bib : 1963ZA.....57...47V .
• Aarseth, Sverre J. (2003). Simulaciones gravitacionales de N cuerpos: herramientas y algoritmos . Cambridge University Press . ISBN 978-0-521-12153-8.
• Bertschinger, Edmund (1998). "Simulaciones de la formación de estructuras en el universo". Annual Review of Astronomy and Astrophysics . 36 (1): 599– 654. Bibcode : 1998ARA & A..36..599B . doi : 10.1146/annurev.astro.36.1.599 .
• Binney, James; Tremaine, Scott (1987). Dinámica galáctica . Princeton University Press . ISBN 978-0-691-08445-9.
• Callahan, Paul B.; Kosaraju, Sambasiva Rao (1992). "Una descomposición de conjuntos de puntos multidimensionales con aplicaciones a k vecinos más cercanos y campos potenciales de n cuerpos (versión preliminar)". STOC '92: Proc. ACM Symp. Theory of Computing . ACM..