¿Cómo "resuelven" los ordenadores el problema de los tres cuerpos?

Question

He investigado un poco y me he enterado de que los ordenadores "resuelven" el problema de los tres cuerpos utilizando "Métodos numéricos para ecuaciones diferenciales ordinarias", pero en realidad no puedo encontrar nada al respecto aparte de Wikipedia. ¿Alguien tiene alguna buena fuente sobre este tema que no sea Wikipedia?

Mis pensamientos:
Actualmente estoy usando simulaciones de tres cuerpos volando uno alrededor del otro, usando la ley gravitacional de Newton, y en un momento aleatorio de la simulación, todo se vuelve caótico. Pensé que esta era la única manera de "resolverlo", pero ¿cómo funciona este método de "Métodos numéricos para ecuaciones diferenciales ordinarias"? ¿Y qué hace realmente el ordenador?

Answer 1

Análisis numérico se utiliza para calcular aproximaciones a cosas: el valor de una función en un punto determinado, dónde está la raíz de una ecuación o las soluciones de un conjunto de ecuaciones diferenciales. Se trata de un tema enorme e importante, ya que en la práctica la mayoría de los problemas reales de matemáticas, ciencia y tecnología no tendrán una solución explícita. solución de forma cerrada (y aunque así fuera, quizá no fuera posible calcular con una precisión infinita; al fin y al cabo, los ordenadores representan los números con una precisión finita). En general, existe un equilibrio entre precisión y velocidad de cálculo.

Para el problema de los tres cuerpos tenemos tres masas puntuales en posiciones iniciales $\mathbf{x}_i(0)$ con velocidades $\mathbf{v}_i(0)$ que queremos calcular para momentos posteriores $t$ . Matemáticamente queremos encontrar la solución del sistema $$\mathbf{x}'_i(t)=\mathbf{v}_i(t),$$ $$\mathbf{v}'_i(t)=\mathbf{f}_i(t)/m_i,$$ $$\mathbf{f}_i(t)=Gm_i \sum_{j\neq i} \frac{m_j(\mathbf{x}_j-\mathbf{x}_i)}{||\mathbf{x}_i-\mathbf{x}_j||^3}.$$

El método obvio es pensar "si avanzamos un pasito $h$ en el tiempo, podemos aproximar todo para que sea lineal", así que hacemos una fórmula donde calculamos el estado en el tiempo $t+h$ del estado en el momento $t$ (y así sucesivamente para $t+2h$ en adelante): $$\mathbf{x}_i(t+h)=\mathbf{x}_i(t)+h\mathbf{v}_i(t),$$ $$\mathbf{v}_i(t+h)=\mathbf{v}_i(t)+h\mathbf{f}_i(t).$$ Esto se llama Método de Euler . Es sencillo pero tiende a ser impreciso; el error por paso es de $\approx O(h^2)$ y tienden a acumularse. Si lo intentas para un problema de dos cuerpos hará que las masas en órbita realicen una órbita de roseta en precesión debido a la acumulación de errores, especialmente cuando se acercan la una a la otra.

Hay una colección de métodos para la resolución numérica de EDO. Se pueden utilizar métodos de orden superior que muestrean las funciones en más puntos y, por tanto, las aproximan mejor. Existen métodos implícitos que en lugar de intentar encontrar un estado en un momento posterior basándose únicamente en el estado actual busque un estado tardío e intermedio autoconsistente. La mayoría de los métodos serios para resolver EDOs también reducirán el tamaño del paso $h$ cuando las fuerzas aumentan durante los encuentros cercanos para garantizar que la precisión siga siendo aceptable. Como ya he dicho, se trata de un gran tema.

Sin embargo, en el caso de las simulaciones mecánicas, es posible que le interesen los métodos diseñados para preservar la energía y otras magnitudes conservadas ( métodos simplécticos - son los que utilizan los profesionales para calcular la órbita a largo plazo). Quizá el más sencillo sea el método de Euler semiimplícito . También existe la Método Verlet y integración leapfrog . Me gusta el método de Euler semiimplícito porque es supersencillo (pero al ser un método de primer orden sigue sin ser terriblemente preciso): $$\mathbf{v}_i(t+h)=\mathbf{v}_i(t)+h\mathbf{f}_i(t),$$ $$\mathbf{x}_i(t+h)=\mathbf{x}_i(t)+h\mathbf{v}_i(t+h).$$ ¿Ves la diferencia? Primero calculas la velocidad actualizada y luego la utilizas para actualizar las posiciones: un pequeño truco, pero de repente las órbitas de 2 cuerpos se comportan bien.

El problema de los tres cuerpos es caótico en un verdadero sentido matemático. Sabemos que hay situaciones en las que pequeñas diferencias en las condiciones iniciales se convierten en diferencias arbitrariamente grandes en las posiciones posteriores (incluso si descartamos los pasos supercercanos entre masas). Por tanto, incluso con una precisión numérica arbitraria, habrá un momento en el que nuestras órbitas calculadas serán totalmente erróneas. El "estilo" general de la trayectoria puede seguir siendo correcto, por lo que está bien jugar con Euler semiimplícito siempre que no se esté planeando ninguna misión espacial basada en los resultados.

Answer 2

La "solución" del problema de los tres cuerpos puede escribirse como un par de ecuaciones diferenciales, \begin{align} \vec{v}&=\frac{\mathrm d\vec{x}}{\mathrm dt}\\ \vec{a}&=\frac{\mathrm d\vec{v}}{\mathrm dt} \end{align} donde esta última suele escribirse en términos de fuerza, $m\vec{a}=\vec{F}$ . A continuación, utilizando la definición de la derivada , $$ \frac{\mathrm df}{\mathrm dx}=\lim_{h\to0}\frac{f(x+h)-f(x)}{h}, $$ las ecuaciones diferenciales con las que empezamos pueden escribirse como las diferencia finita ecuaciones , \begin{align} \vec x(t+\mathrm dt)&\simeq \vec x(t)+\vec v\,\mathrm dt \\ \vec v(t+\mathrm dt)&\simeq \vec v(t)+\vec a\,\mathrm dt \end{align} lo que se hace suponiendo que $\mathrm dt$ es "sólo pequeño" y no infinitesimal.

Es estos que un ordenador utiliza para "resolver" el problema de los tres cuerpos: dada una condición inicial para cada uno de los cuerpos, el ordenador avanza iterativamente en el tiempo utilizando las ecuaciones de diferencias finitas anteriores, siendo la fuerza la variable $n$ -fuerza del cuerpo .

Como se menciona en los comentarios, el modelo simplista anterior se denomina Integración de Euler que no se adapta en absoluto a este problema porque no conserva la energía. Una opción mejor es Velocidad Verlet que se describe en otros posts en physics.SE (Recomiendo este post mío ya que ofrece algunos detalles decentes, aunque breves, de la aplicación).

Un documento razonablemente accesible sobre el problema de los tres cuerpos, que incluye una discusión sobre los algoritmos numéricos, es Musielak y Quarles (2014) . Para un análisis más amplio de las técnicas de orden superior más útiles, como se sugiere por homocomputeris Ver Farrés y otros (2012)

Answer 3

Resulta que conociendo las trayectorias por pares se puede determinar la ecuación de movimiento de los cuerpos. Las ecuaciones de movimiento se escriben como $$\frac{d^2 \vec r^k}{ds^2}=-G\sum_{n=1,n\ne k}^N \frac{m_n(\vec r^k-\vec r^n)}{|\vec r^k-\vec r^n|^3}\tag{1}$$ Resolvemos el problema auxiliar de la interacción de pares de cuerpos con masa inercial reducida $$\frac{m_n m_k}{\sum_{k=1}^N m_k} \frac{d^2 \vec R^{kn}}{ds^2}=-G \frac{m_n m_k\vec R^{kn}}{|\vec R^{kn}|^3 }\tag{2}$$ Estas ecuaciones difieren en las distintas condiciones iniciales. Restando de la ecuación (2) la ecuación (1), obtenemos $$\frac{d^2 R_0^k}{ds^2}-\frac{d^2 r^k}{ds^2}= G\sum_{n=1,n\ne k}^N \frac{m_n(\vec r^k-\vec r^n)}{|\vec r^k-\vec r^n|^3}- G\sum_{n=1,n\ne k}^N \frac{m_n\vec R^{kn}}{|\vec R^{kn}|^3}=0,\vec R^{kn}=\vec r^k-\vec r^n$$ Dónde $\vec R_0^k=\sum_{n=1 n \ne k}^N \frac{m_n \vec R^{kn}}{\sum_{p=1}^N m_p}$ Obtener la fórmula $$\frac{d^2 \vec R_0^k}{ds^2}=\frac{d^2 \vec r^k}{ds^2}$$ Integrando esta igualdad, obtenemos la ecuación de movimiento para cada uno de los N cuerpos. $$\vec r_0^k=\frac{d\vec r^k(0)}{ds}s+\vec r^k(0)+\vec R_0^k(s)- \frac{d\vec R^k(0)}{ds}s-\vec R^k(0)= \vec R_0^k(s)$$ El movimiento de cada cuerpo está determinado por el movimiento del centro de inercia del sistema de cuerpos pares. Las condiciones iniciales se determinan a partir de la condición $$\frac{d\vec r^k(0)}{ds}s+\vec r^k(0)= \frac{d\vec R^k(0)}{ds}s+\vec R^k(0)\tag{3}$$ Seleccionamos el sistema de coordenadas $\sum_{n=1}^N m_n \vec r^n(0)=0,\sum_{n=1}^N m_n \frac{d \vec r^n(0)}{ds}=0$ y entonces la condición (3) se cumple idénticamente. Obtendremos $ \vec R^{kn}(0)=\vec r^k(0)- \vec r^n(0), \frac{d \vec R^{kn}}{ds}(0)=\frac{d \vec r^k}{ds}(0)-\frac{d \vec r^n}{ds}(0)$ La energía de la trayectoria par del sistema se calcula mediante la fórmula $$E_{kn}=\frac{m_k m_n}{\sum_{q=1}^N m_q} (\frac{d \vec R^{kn}}{ds}(0))^2/2-\frac{G m_k m_n}{|\vec R^{kn}(0)|}+\frac{M_{kn}^2\sum_{q=1}^N m_q}{2 m_k m_n |R^{kn}|^2}=\frac{m_k m_n}{\sum_{q=1}^N m_q} (\frac{d \vec R^{kn}}{ds}(0))^2/2-\frac{G m_k m_n}{|\vec R^{kn}(0)|}+\frac{m_k m_n}{2 \sum_{q=1}^N m_q |R^{kn}(0)|^2}(\frac{d \vec R^{kn}}{ds}(0)\times \vec R^{kn}(0))^2$$ El momento del impulso viene determinado por la fórmula producto vectorial $$M_{kn}=\frac{m_k m_n}{\sum_{q=1}^N m_q} \frac{d \vec R^{kn}}{ds}(0)\times \vec R^{kn}(0)$$ En este caso, la solución construida se sustituye en las trayectorias calculadas $\vec r_0^k-\vec r_0^p $ de los planetas y utilizó $\vec R^{kp}=r^k-r^p$ obtenemos $\vec R^{kp}=\vec r_0^k-\vec r_0^p$ . De hecho $$ r_0^k-r_0^p =\sum_{n=1}^N \frac{m_n \vec R^{kn}}{\sum_{q=1}^N m_q}-\sum_{n=1}^N \frac{m_n \vec R^{pn}}{\sum_{q=1}^N m_q}=\sum_{n=1}^N \frac{m_n (\vec r^k-\vec r^n)}{\sum_{q=1}^N m_q}-\sum_{n=1}^N \frac{m_n (\vec r^p-\vec r^n)}{\sum_{q=1}^N m_q}=\vec r^k-\vec r^p=R_{kp}$$ Fórmula simplificada $$\vec r_0^u=\sum_{n=1}^N \frac{m_n \vec R^{un}(s)}{\sum_{q=1}^N m_q |R^{un}(s)|} \frac{p_{un}}{1+e_{un}cos(\varphi_{un}-\varphi_{un0})}$$ El ángulo se calcula mediante la fórmula $\varphi_{un}-\varphi_{un0}=\int_0^s \frac{M_{un}}{|\vec R^{pn}(u)|^2}du$ La magnitud $$\frac{R^{un}(\varphi)}{|R^{un}(\varphi)|}=\frac{\vec R^{un}(0)}{|\vec R^{un}(0)|}cos\varphi+\frac{\vec R^{un}(0)\times \vec M_{un}}{|\vec R^{un}(0)\times \vec M_{un}|}sin\varphi $$ del ángulo $\varphi$ común a todas las trayectorias pares corresponde a la dirección $\frac{R^{un}(\varphi)}{|R^{un}(\varphi)|}$ con ortos individuales $$\frac{\vec R^{un}(0)}{|\vec R^{un}(0)|},\frac{\vec R^{un}(0)\times \vec M_{un}}{|\vec R^{un}(0)\times \vec M_{un}|} $$ . Tenemos la fórmula final para la trayectoria de los cuerpos $$\vec r_0^u(\varphi)=\sum_{n=1}^N \frac{m_n}{\sum_{q=1}^N m_q}[\frac{\vec R^{un}(0)}{|\vec R^{un}(0)|}cos\varphi+\frac{\vec R^{un}(0)\times \vec M_{un}}{|\vec R^{un}(0)\times \vec M_{un}|}sin\varphi]\frac{p_{un}}{1+e_{un}cos\varphi}$$ Pero a falta de una solución infinita aunque sea por poco tiempo, es necesario suponer $$|\vec R^{un}(\varphi)|=\frac{p_{un}}{1+e_{un}cos\varphi}>0,e_{un}<1$$ Esto impone una restricción al valor de los parámetros; la energía de la interacción de pares debe ser negativa. Así, se excluyen los cuerpos de masa pequeña, cuya contribución a la energía del sistema no es significativa. De acuerdo con la fórmula de la interacción de pares, el radio varía según la ley véase Landau Lifshits volumen 1 $$r=a(ecosh\xi-1);t=\sqrt{\frac{a^3}{G\sum_{q=1}^N m_q}}(esinh\xi-\xi),r\approx c\sqrt{\frac{r_g}{a}}t$$ Tras sencillas transformaciones, esta fórmula se reduce a $$r=\sqrt{V^2(0)+\frac{(\vec V(0)\times \vec R(0))^2}{R(0)^2}-\frac{2G\sum_{q=1}^N m_q}{R(0)}}t$$ Este aumento lineal del radio desplaza el centro de gravedad del sistema. Si una partícula se desplaza hacia el infinito, todo el sistema es arrastrado y el centro de gravedad se desplaza. La segunda velocidad cósmica corresponde a un valor positivo del radicando, o energía positiva del sistema.

Answer 4

Hace unos años monté una maqueta del sistema solar interior: el Sol, la Luna terrestre y los planetas desde Mercurio hasta Marte.

Utilicé una variación Leapfrog de Rk4 y obtuve resultados decentes, períodos orbitales coherentes con las configuraciones iniciales que se mantuvieron estables en el tiempo. Incluso algunas de las cosas malas que obtuve creo que sugerían buenas prácticas numéricas.

Inicializaría mis planetas con momento inicial, con el sol teniendo momento inicial cero. El impulso de los planetas no se anulaba, de modo que se producía una extraña deriva neta en una dirección aproximadamente rectilínea que tardé un tiempo en comprender.

Además de solucionar eso, una característica que me pareció importante fue el orden de aplicación de los cambios de posición y velocidad.

No querrás iterar a través de tus cuerpos, calcular la fuerza neta en una pasada cuerpo por cuerpo, y luego aplicar tus cambios dinámicos cuerpo por cuerpo. De lo contrario, por ejemplo, actualizar la posición de Venus, manteniendo la posición de Marte sin cambios, mientras que usted está actualizando la posición de la Tierra. Tus pasos temporales se desincronizan.

Almacene en algún lugar de la memoria los nuevos desplazamientos y velocidades a medida que itera a través de los cuerpos, y sólo después de encontrar su dinámica en el tick actual, ejecute y aplique sus cambios. Mis disculpas y enhorabuena si te diste cuenta en tu primer intento. Me hizo tropezar un poco.

Answer 5

Una forma de simular tres cuerpos es utilizando incrementos de paso de tiempo, en los que se pretende que cada cuerpo acelera a un ritmo constante durante un breve periodo de tiempo. Se comienza con los vectores iniciales de posición, velocidad y masas de los cuerpos, y luego se utilizan para calcular el vector de aceleración de cada cuerpo. En Física Newtoniana el vector aceleración gravitatoria del Cuerpo A respecto del Cuerpo B viene dado por la ecuación $$\vec{a_{GA}}=-\frac{G\left(\vec{A}-\vec{B}\right)M_B}{||\vec{A}-\vec{B}||^3}$$ con $G$ que es la constante gravitatoria, $\vec{A}$ siendo el vector de posición del Cuerpo A, $\vec{B}$ que es el vector de posición del cuerpo B, $M_B$ siendo la masa del cuerpo B, y $\vec{a_{GA}}$ siendo el Vector de Aceleración Gravitacional del Cuerpo A. Para encontrar el vector de aceleración gravitacional total del Cuerpo A se suman cada uno de los vectores de aceleración gravitacional del cuerpo A de cada uno de los otros cuerpos.

Puede encontrar los vectores de posición y los vectores de velocidad al final del paso de tiempo $n$ utilizando las fórmulas generales $$\vec{x_n}=\frac{1}{2}\vec{a_{n-1}}s^2+\vec{v_{n-1}}s+\vec{x_{n-1}}$$ $$\vec{v_n}=\vec{a_{n-1}}s+\vec{v_{n-1}}$$ con $\vec{a_{n-1}}$ siendo el vector de aceleración para el paso de tiempo $n-1$ , $\vec{v_{n-1}}$ que es el vector de velocidad para el paso de tiempo $n-1$ , $\vec{x_{n-1}}$ siendo el vector de posición $n-1$ , $s$ siendo la longitud de los incrementos de los pasos temporales, $\vec{x_n}$ siendo el vector de posición para el paso de tiempo $n$ y $\vec{v_n}$ siendo el vector de velocidad para el paso de tiempo $n$ . En cada paso temporal también se recalculan los vectores de aceleración, ya que no son constantes.

Este método de simulación $n$ funciona para cualquier ley de fuerza, lo que significa que se puede utilizar para simular $n$ organismos para una ley de la fuerza $F=f(r)$ con $f(r)$ siendo alguna función de la distancia, aunque en el caso de la ley de fuerza $F=f(r)$ la ecuación que especifiqué anteriormente para calcular la aceleración del cuerpo A desde el cuerpo B se convierte en $$\vec{a_{GA}}=-\frac{\vec{A}-\vec{B}}{||\vec{A}-\vec{B}||}GM_Bf\left(||\vec{A}-\vec{B}||\right)$$ pero las demás fórmulas que he mencionado, incluida la parte de sumar los vectores de aceleración, siguen siendo las mismas.