Simulación de n cuerpos

Una simulación de n cuerpos de la formación cosmológica de un cúmulo de galaxias en un universo en expansión.

En física y astronomía, una simulación de n cuerpos es una simulación de un sistema dinámico de partículas, generalmente bajo la influencia de fuerzas físicas, como la gravedad (ver problema de los n cuerpos para otras aplicaciones). Las simulaciones de n cuerpos son herramientas ampliamente utilizadas en astrofísica, desde la investigación de la dinámica de sistemas de pocos cuerpos como el sistema Tierra-Luna-Sol hasta la comprensión de la evolución de la estructura a gran escala del universo.​ En cosmología física, las simulaciones de n cuerpos se utilizan para estudiar procesos de formación de estructuras no lineales, como filamentos y halos de galaxias, debido a la influencia de la materia oscura. También se utilizan simulaciones directas de n cuerpos para estudiar la evolución dinámica de los cúmulos estelares.

Naturaleza de las partículas

Las 'partículas' tratadas por la simulación pueden corresponder o no a objetos físicos de naturaleza definida por partículas. Por ejemplo, una simulación de n cuerpos de un cúmulo de estrellas podría estar compuesto de n partículas, por lo que cada partícula representa una estrella y tiene un significado físico. Por otro lado, en el caso de una simulación de una nube de gas, resulta imposible representar cada átomo o molécula de gas por una partícula, ya que esto requeriría del orden de 10²³ partículas por cada mol de material (ver constante de Avogadro). Debido a limitaciones de la computación actual, una sola 'partícula' debe representar una cantidad mucho mayor de gas (a menudo implementado utilizando Hidrodinámica de Partículas Suavizadas). Esta cantidad no necesita tener ningún significado físico, sino que debe elegirse como un compromiso entre precisión y requisitos informáticos manejables.

Simulación de Materia Oscura

La materia oscura toma un papel importante en la formación de galaxias. La ecuación de Boltzmann sin colisiones describe la evolución de la densidad f de las partículas de materia oscura respecto al tiempo (en el espacio de fases):

$${\frac {df}{dt}}={\frac {\partial f}{\partial t}}+\mathbf {v} \cdot \nabla f-{\frac {\partial f}{\partial \mathbf {v} }}\cdot \nabla \Phi$$

En esta ecuación, $\mathbf {v}$ es la velocidad y Φ es el potencial gravitacional obtenido por la ecuación de Poisson. Estas dos ecuaciones juntas se resuelven en un universo que se rige por las ecuaciones de Friedmann tras determinar las condiciones iniciales de las partículas de materia oscura. El método convencional para inicializar posiciones y velocidades implica mover las partículas dentro de una red cartesiana uniforme o una configuración de partículas similar al vidrio.​ Esto se puede conseguir mediante una aproximación de teoría lineal o una teoría de perturbaciones de bajo orden.​

Simulaciones de n cuerpos gravitacionales directos

Simulación de n cuerpos de 400 objetos con parámetros similares a los de los planetas del Sistema Solar .

En las simulaciones gravitacionales directas de n cuerpos, estos se encuentran bajo la mutua influencia de sus fuerzas gravitacionales. Las respectivas ecuaciones se integran numéricamente sin aproximaciones ni simplificaciones. Estos cálculos se utilizan en casos donde las interacciones entre objetos individuales, como estrellas o planetas, son importantes para la evolución del sistema.

Las primeras simulaciones gravitacionales directas de n cuerpos fueron llevadas a cabo por Erik Holmberg en el Observatorio de Lund en 1941, determinando las fuerzas entre estrellas en galaxias influidas mutuamente por su gravitación mediante la equivalencia matemática entre la propagación de la luz y la interacción gravitacional. Colocando bombillas en las posiciones de las estrellas y midiendo los flujos de luz direccionales con una célula fotoeléctrica, las ecuaciones de movimiento se pueden integrar con una complejidad de O(N).​ Las primeras simulaciones puramente por cálculos fueron realizadas por Sebastian von Hoerner en el Astronomisches Rechen-Institut de Heidelberg, Alemania. Sverre Aarseth de la Universidad de Cambridge (Reino Unido) ha dedicado toda su vida científica al desarrollo de una serie de códigos de n cuerpos altamente eficientes para aplicaciones astrofísicas que utilizan pasos de tiempo adaptativos (jerárquicos), un esquema vecino de Ahmad-Cohen y la regularización de encuentros cercanos. La regularización es un truco matemático para eliminar la singularidad en la ley de gravitación newtoniana de dos partículas que se acercan arbitrariamente entre sí. Los códigos de Sverre Aarseth se utilizan para estudiar la dinámica de cúmulos estelares, sistemas planetarios y núcleos galácticos.^{[cita requerida]}

Simulaciones de relatividad general

Muchas simulaciones son lo suficientemente grandes como para que los efectos de la relatividad general sean significativos al establecer una cosmología de Friedmann-Lemaitre-Robertson-Walker. Esto se incorpora en la simulación como una medida evolutiva de la distancia (o factor de escala) en un sistema de coordenadas en movimiento, lo que hace que las partículas se desaceleren en las coordenadas en movimiento (así como debido al desplazamiento al rojo de su energía física). Sin embargo, las contribuciones de la relatividad general y la velocidad finita de la gravedad pueden ignorarse, ya que las escalas de tiempo dinámicas típicas son largas en comparación con el tiempo de cruce de la luz para la simulación, y la curvatura espacio-temporal inducida por las partículas y las velocidades de las partículas son pequeñas. Las condiciones de contorno de estas simulaciones cosmológicas suelen ser periódicas (o toroidales), de modo que un borde del volumen de simulación coincide con el borde opuesto.

Optimizaciones de cálculo

Las simulaciones de N cuerpos son simples en principio, porque implican simplemente integrar las 6 N ecuaciones diferenciales ordinarias que definen los movimientos de las partículas en la gravedad newtoniana . En la práctica, el número N de partículas involucradas suele ser muy grande (las simulaciones típicas incluyen muchos millones, la simulación Millennium incluyó diez mil millones) y el número de interacciones entre partículas que necesitan calcularse aumenta del orden de N², y por lo tanto, La integración de las ecuaciones diferenciales puede resultar costosa o incluso imposible desde el punto de vista computacional. Por lo tanto, comúnmente se utilizan una serie de mejoras.

La integración numérica generalmente se realiza en pequeños intervalos de tiempo utilizando un método como la integración por salto . Sin embargo, toda integración numérica conduce a errores. Los pasos más pequeños dan menos errores pero se ejecutan más lentamente. La integración de Leapfrog es aproximadamente de segundo orden en el paso de tiempo, otros integradores, como los métodos de Runge-Kutta, pueden tener una precisión de cuarto orden o mucho mayor.

Uno de los refinamientos más simples es que cada partícula lleva consigo su propia variable de paso de tiempo, de modo que las partículas con tiempos dinámicos muy diferentes no tienen que evolucionar hacia adelante al ritmo del tiempo más corto.

Existen dos esquemas de aproximación básicos para disminuir el tiempo de cálculo de tales simulaciones. Estos pueden reducir la complejidad computacional a O (N log N) o mejor, con pérdida de precisión.

Métodos de árbol

En los métodos de árbol, como la simulación de Barnes-Hut, generalmente se usa un octree para dividir el volumen en celdas cúbicas y solo las interacciones entre partículas de celdas cercanas deben tratarse individualmente; Las partículas en células distantes pueden tratarse colectivamente como una única partícula grande centrada en el centro de masa de la célula distante (o como una expansión multipolar de bajo orden). Esto puede reducir drásticamente el número de interacciones de pares de partículas que deben calcularse. Para evitar que la simulación se vea saturada por el cálculo de las interacciones entre partículas, las celdas deben refinarse a celdas más pequeñas en 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 O óptimo ( n registro n ) tiempo por iteración con dimensión fija.

Método de malla de partículas

Otra posibilidad es el método de malla de partículas en el que el espacio se discretiza en una malla y, a los efectos de calcular el potencial gravitacional, se supone que las partículas están divididas entre los vértices circundantes de 2x2 de la malla. Encontrar la energía potencial Φ es fácil, porque la ecuación de Poisson

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

donde G es la constante de Newton y $\rho$ es la densidad (número de partículas en los puntos de la malla), es trivial de resolver usando la transformada rápida de Fourier para ir al dominio de la frecuencia donde la ecuación de Poisson tiene la forma simple

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

dónde ${\vec {k}}$ es el número de onda comóvil y los sombreros denotan transformadas de Fourier. Desde ${\vec {g}}=-{\vec {\nabla }}\Phi$, el campo gravitacional ahora se puede encontrar multiplicando por $-i{\vec {k}}$ y calcular la transformada inversa de Fourier (o calcular la transformada inversa y luego utilizar 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 o alguna otra técnica (como la combinación con un árbol o un algoritmo simple partícula-partícula) 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.

Optimizaciones de casos especiales

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

La gente suele decidir poner un satélite en una órbita congelada . La trayectoria de un satélite que orbita cerca de la Tierra se puede modelar con precisión a partir de la órbita elíptica de dos cuerpos alrededor del centro de la Tierra y agregando pequeñas correcciones debido al achatamiento de la Tierra, la atracción gravitacional del Sol y la Luna y la resistencia atmosférica., etc. Es posible encontrar una órbita congelada sin calcular la trayectoria real del satélite.

La trayectoria de un pequeño planeta, cometa o 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 agregando pequeñas correcciones de la atracción gravitacional 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 se pueden calcular directamente. No es necesario calcular la trayectoria real de ninguna partícula en particular como paso intermedio. Tales características incluyen la estabilidad de Lyapunov, el tiempo de Lyapunov, varias medidas de la teoría ergódica, etc.

Sistemas de dos partículas

Aunque hay millones o miles de millones de partículas en las simulaciones típicas, normalmente corresponden a una partícula real con una masa muy grande, normalmente 10⁹ masas solares. Esto puede introducir problemas con interacciones de corto alcance entre las partículas, como la formación de sistemas binarios de dos partículas. Como se supone que las partículas representan una gran cantidad de partículas de materia oscura o grupos de estrellas, estos binarios no son físicos. Para evitar esto, se utiliza una ley de fuerza newtoniana suavizada, que no diverge como el radio del cuadrado inverso en distancias cortas. La mayoría de las simulaciones implementan esto de forma bastante natural ejecutando las simulaciones en celdas de tamaño finito. Es importante implementar el procedimiento de discretización de tal manera que las partículas siempre ejerzan una fuerza evanescente sobre sí mismas.

Suavizamiento

El suavizamiento 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 llega al infinito). Esto se obtiene modificando el potencial gravitacional regularizado de cada partícula como

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

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

Resultados de simulaciones de n cuerpos

Las simulaciones de n cuerpos arrojan resultados sobre la distribución de la materia oscura a gran escala y la estructura de los halos de materia oscura. Según simulaciones de materia oscura fría, la distribución general de la materia oscura a gran escala no es del todo uniforme. En cambio, muestra una estructura que se asemeja a una red, que consta de vacíos, paredes, filamentos y halos. Además, las simulaciones muestran que la relación entre la concentración de los halos y factores como la masa, el espectro de fluctuación inicial y los parámetros cosmológicos está relacionada con el tiempo real de formación de los halos.​ En particular, los halos con menor masa tienden a formarse antes y, como resultado, tienen concentraciones más altas debido a la mayor densidad del Universo en el momento de su formación. Se ha descubierto que las formas de los halos no son perfectamente esféricas. Por lo general, los halos se alargan y se prolongan cada vez más hacia sus centros. Sin embargo, las interacciones entre la materia oscura y los bariones afectarían la estructura interna de los halos de materia oscura. Se necesitan simulaciones que modelen tanto la materia oscura como los bariones para estudiar estructuras a pequeña escala.

Incorporando bariones, leptones y fotones en simulaciones.

Muchas simulaciones simulan únicamente materia oscura fría y, por tanto, incluyen únicamente la fuerza gravitacional. La incorporación de bariones, leptones y fotones en las simulaciones aumenta drásticamente su complejidad y, a menudo, es necesario realizar simplificaciones radicales de la física subyacente. Sin embargo, esta es un área extremadamente importante y muchas simulaciones modernas ahora están tratando de comprender los procesos que ocurren durante la formación de galaxias y que podrían explicar el sesgo de las galaxias.

Complejidad computacional

Reif y Tate​ demuestran que si el problema de alcanzabilidad de n cuerpos se define de la siguiente manera: dados n cuerpos que satisfacen una ley de potencial electrostático fija, se determina si un cuerpo alcanza una bola de destino en un período de tiempo determinado donde requerimos un poli(n) bits de precisión y el tiempo objetivo es poli(n) está en PSPACE.

Por otro lado, si la pregunta es si el cuerpo finalmente llega a la bola de destino, el problema es difícil en PSPACE. Estos límites se basan en límites de complejidad similares obtenidos para el trazado de rayos.

Simulaciones de ejemplo

Código repetitivo común

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; ingenuo al implicar que las únicas fuerzas que actúan sobre los cuerpos en órbita es la fuerza gravitacional que ejercen entre sí. En lenguajes de programación orientados a objetos, como C++, algún 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 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 sería por atracción o repulsión por interacción de campos eléctricos. De todos modos, la aceleración de una partícula es el resultado de la suma de los vectores de fuerza, divididos por la masa de la partícula:

$${\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, que en código optimizado permite una fácil asignación de memoria y predicción de los recursos consumidos; como se ve en el siguiente código C++:

struct Vector3
{
    double e[3] = { 0 };

    Vector3() {}
    ~Vector3() {}

    inline Vector3(double e0, double e1, double e2)
    {
        this->e[0] = e0;
        this->e[1] = e1;
        this->e[2] = e2;
    }
};

struct OrbitalEntity
{
    double e[7] = { 0 };

    OrbitalEntity() {}
    ~OrbitalEntity() {}

    inline OrbitalEntity(double e0, double e1, double e2, double e3, double e4, double e5, double e6)
    {
        this->e[0] = e0;
        this->e[1] = e1;
        this->e[2] = e2;
        this->e[3] = e3;
        this->e[4] = e4;
        this->e[5] = e5;
        this->e[6] = e6;
    }
};

Tenga en cuenta que OrbitalEntity contiene 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 $\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 $\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 $\left[0\;0\;1\right]$
• ${\textstyle e_{3}={\dot {x}}}$, la proyección del vector velocidad de los objetos en el espacio cartesiano a lo largo $\left[1\;0\;0\right]$
• ${\textstyle e_{4}={\dot {y}}}$, la proyección del vector velocidad de los objetos en el espacio cartesiano a lo largo $\left[0\;1\;0\right]$
• ${\textstyle e_{5}={\dot {z}}}$, la proyección del vector velocidad de los objetos en el espacio cartesiano a lo largo $\left[0\;0\;1\right]$

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

Inicialización de parámetros de simulación.

Comúnmente, las simulaciones de N cuerpos serán sistemas basados en algún tipo de ecuaciones de movimiento ; De estos, la mayoría dependerá de alguna configuración inicial para "iniciar" la simulación. En sistemas como los que dependen de algún potencial eléctrico o gravitacional, la fuerza sobre una entidad de simulación es independiente de su velocidad. Por lo tanto, para sembrar las fuerzas de la simulación, simplemente se necesitan posiciones iniciales, pero esto no permitirá la propagación: se requieren velocidades iniciales. Considere un planeta que orbita alrededor de una estrella: no tiene movimiento, pero está sujeto a la atracción gravitacional hacia 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, $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 $t_{n+1}$, el cambio de posición resultante es significativamente diferente debido a la dependencia inherente de la propagación con la velocidad. En mecanismos de propagación básicos, como el método simpléctico de Euler que se utilizará a continuación, la posición de un objeto en $t_{n+1}$ sólo depende de su velocidad en $t_{n}$. 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 lograr una simulación similar a la del sistema solar tomando distancias promedio de masas puntuales equivalentes a planetas desde una estrella central. Para mantener el código simple, se utiliza un enfoque no riguroso basado en semiejes mayores y velocidades medias. Se debe reservar espacio de memoria para estos cuerpos antes de configurarlos; Para permitir la escalabilidad, se puede utilizar un 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 };   // a star similar to the sun
orbital_entities[1] = { 57.909e9,0.0,0.0,   0.0,47.36e3,0.0,  0.33011e24 }; // a planet similar to mercury
orbital_entities[2] = { 108.209e9,0.0,0.0,  0.0,35.02e3,0.0,  4.8675e24 };  // a planet similar to venus
orbital_entities[3] = { 149.596e9,0.0,0.0,  0.0,29.78e3,0.0,  5.9724e24 };  // a planet similar to earth
orbital_entities[4] = { 227.923e9,0.0,0.0,  0.0,24.07e3,0.0,  0.64171e24 }; // a planet similar to mars
orbital_entities[5] = { 778.570e9,0.0,0.0,  0.0,13e3,0.0,     1898.19e24 }; // a planet similar to jupiter
orbital_entities[6] = { 1433.529e9,0.0,0.0, 0.0,9.68e3,0.0,   568.34e24 };  // a planet similar to saturn
orbital_entities[7] = { 2872.463e9,0.0,0.0, 0.0,6.80e3,0.0,   86.813e24 };  // a planet similar to uranus
orbital_entities[8] = { 4495.060e9,0.0,0.0, 0.0,5.43e3,0.0,   102.413e24 }; // a planet similar to neptune

donde N_ASTEROIDS es 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 simulaciones es establecer los rangos de tiempo de la simulación, $t_{0}$ a $t_{\text{end}}$, así como el paso de tiempo incremental $dt$ que hará avanzar la simulación:

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

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

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

Propagación

Una simulación completa puede constar de cientos, miles, millones, miles de millones o, a veces, billones de pasos de tiempo. En el nivel elemental, cada paso de tiempo (para simulaciones con partículas que se mueven debido a fuerzas ejercidas sobre ellas) implica

• calcular las fuerzas sobre cada cuerpo
• calculando las aceleraciones de cada cuerpo ( ${\vec {a}}$ )
• calculando las velocidades de cada cuerpo ( ${\vec {v}}_{n}={\vec {v}}_{n-1}+{\vec {a}}_{n}$ )
• calculando la nueva posición de cada cuerpo ( ${\vec {r}}_{n+1}={\vec {r}}_{n}+{\vec {v}}_{n}$ )

Lo anterior se puede implementar de manera bastante simple con un bucle while que continúa mientras $t$ existe en el rango antes mencionado:

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

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

Bibliografía

• von Hoerner, Sebastian (1960). «Die numerische Integration des n-Körper-Problemes für Sternhaufen. I». Zeitschrift für Astrophysik (en alemán) 50: 184. Bibcode:1960ZA.....50..184V. 
• von Hoerner, Sebastian (1963). «Die numerische Integration des n-Körper-Problemes für Sternhaufen. II». Zeitschrift für Astrophysik (en alemán) 57: 47. Bibcode:1963ZA.....57...47V. 
• Aarseth, Sverre J. (2003). Gravitational N-body Simulations: Tools and Algorithms. Cambridge University Press. ISBN 978-0-521-12153-8. 
• Bertschinger, Edmund (1998). «Simulations of structure formation in the universe». 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). Galactic Dynamics. Princeton University Press. ISBN 978-0-691-08445-9. 
• Callahan, Paul B.; Kosaraju, Sambasiva Rao (1992). «A decomposition of multidimensional point sets with applications to k-nearest-neighbors and n-body potential fields (preliminary version)». STOC '92: Proc. ACM Symp. Theory of Computing. ACM. .

• Datos: Q85787439