Numéricos: ¿Cómo renormalizo el siguiente ODE?

Esta pregunta es más acerca de cómo abordar un problema numéricamente.

En un pequeño proyecto quería simular el movimiento coorbital de Janus y Epimetheus. Esto es básicamente un problema de tres cuerpos. Elijo que Saturno se fije en el origen, dejemos que $r_1$ y $r_2$ sean los vectores de ubicación de janus y epimeteo, respectivamente. Dado que el efecto ocurre cuando Janus y Epimeteo están muy juntos, elegí coordenadas relativas para una mejor resolución, es decir, $r=r_1-r_2$ y $R=r_1+r_2$ . Ahora obtengo las siguientes ecuaciones de movimiento:

\frac{d^{2}}{d t^{2}} (\binom{R}{r}) = - G (m_{2} \pm m_{1}) \frac{R}{R^{3}} - 4 M G (\frac{r + R}{(r + R)^{3}} \mp \frac{r - R}{(r - R)^{3}})

$\frac {d^2}{dt^2} \binom{R}{r} = - G (m_2\pm m_1) \frac R {R^3} - 4 M G \left(\frac {r+R}{(r+R)^3} \mp \frac {r-R}{(r-R)^3}\right )$

donde corresponde a las masas de las lunas, es la masa de Saturno y la constante gravitacional. El problema surge cuando trato de resolver esto numéricamente. Uno tiene que lidiar con valores de magnitudes completamente diferentes, es decir, y . Y , se encuentran en las regiones de 0 a 150.000. $m_i$ $M$ $G$ $M \sim e^{28}$ $m_i \sim e^{17}$ $r$ $R$

Para ser honesto, no estoy seguro si este es el foro para discutir tales problemas numéricos.

Más información:

El código está escrito en Matlab y utilizo un solucionador ODE estándar para obtener el resultado. Sin embargo, esto se está desmoronando porque el tamaño del paso no se puede reducir con la precisión de la máquina. (Creo que esto no es sorprendente porque uno tiene que lidiar con los órdenes de magnitud ya mencionados).

ode GertVdE
fuente

¿Estás ejecutando esta simulación en unidades SI? Como mínimo, debe dividir todo por algún factor de , de modo que pueda eliminar algunos órdenes de magnitud.

G * m_{2}

$G*m_{2}$

Hola, esto, pero todavía no funciona ... Se producen los mismos problemas que antes. :(

Debe establecer su unidad de masa en una de las masas de la luna, y sus unidades de longitud / tiempo para establecer las cosas en 1. Nada debe ser menor que 1/100 si lo escribe bien. No hay necesidad de un solucionador de venta libre. Escriba el código para hacerlo usted mismo, donde controle el tamaño de los pasos. Las interrupciones en el tamaño de pasos con este tipo de potenciales pueden ocurrir en colisiones, donde el solucionador intentará reducir el tamaño de pasos hasta la convergencia, y en la colisión no hay convergencia. Debe asegurarse de que las órbitas no sean colineales, por lo que debe ver la simulación. No puede obtener una respuesta como es.

Ron Maimon

Evite las abreviaturas en el título. DGL = Diferencial de ángulo?

¿Qué solucionador de ODE estándar está utilizando?

Geoff Oxberry

Respuestas:

Su enfoque actual arruina la estabilidad numérica; de hecho, probablemente pierdas resolución de esta manera.

Tome como coordenadas para cada satélite sus variables Kepler y el ángulo del plano que contiene la posición del satélite, la velocidad y el origen. Las ecuaciones diferenciales en ausencia de interacción entre los satélites son trivialmente simples, y solo la interacción se vuelve algo complicada. Como la interacción es pequeña si los satélites están muy lejos, la dinámica resultante debería ser numéricamente estable.

Arnold Neumaier
fuente

En lugar de utilizar un solucionador de ODE "clásico" (rígido), puede utilizar algoritmos dedicados para la integración numérica geométrica. Consulte, por ejemplo, este libro y los códigos GNI que puede encontrar en el sitio web de Ernst Hairer .

GertVdE
fuente

¿Qué tal si tiene tres pasos en su simulación:

actualizar la posición de Janus calculando la fuerza de Janus - Saturno.
actualizar la posición de Epimeteo calculando la fuerza de Epimeteo - Saturno.
actualizar la posición de Janus y Epimetheus calculando la fuerza de Janus - Epimetheus.

Posiblemente utilizando mejores tiempos para # 3.

No estoy seguro si esto ayudará. Supongo que el verdadero problema es que la magnitud de la fuerza es diferente en la luna - luna y la luna - caso de Saturno, excepto si las lunas están cerca.

Alternativamente:

si las lunas se cierran, calcule unas lunas aproximadas: la fuerza de Saturno usando su centro de masa vectorial y actualice ambas posiciones con el mismo vector.
si están muy separados, actualícelos por separado.
como antes.

¡La mejor de las suertes!

fuente