Numéricos: Como renormalizo o seguinte ODE

Esta pergunta é mais sobre como resolver um problema numericamente.

Em um pequeno projeto, eu queria simular o movimento coorbital de Janus e Epimetheus. Este é basicamente um problema de três corpos. Eu escolhi Saturno para ser corrigido na origem, sejam e os vetores de localização de janus e epimetheus, respectivamente. Como o efeito ocorre quando Janus e Epimetheus estão muito próximos, escolhi coordenadas relativas para uma melhor resolução, ou seja, e . Agora eu recebo as seguintes equações de movimento: $r_1$ $r_2$ $r=r_1-r_2$ $R=r_1+r_2$

\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 )$

onde $m_i$ corresponde às massas das luas, $M$ é a massa de Saturno e $G$ a constante gravitacional. O problema surge quando tento resolver isso numericamente. É preciso lidar com valores de magnitudes completamente diferentes, isto é, $M \sim e^{28}$ e $m_i \sim e^{17}$ . E $r$ , $R$ estão nas regiões de 0 a 150.000.

Para ser sincero, não tenho certeza se este é o fórum para discutir esses problemas numéricos.

Mais Informações:

O código está escrito no Matlab e eu uso um solucionador de ODE padrão para obter o resultado. No entanto, isso está quebrando porque o tamanho da etapa não pode ser reduzido com precisão da máquina. (Acho que isso não é surpreendente, porque é preciso lidar com as ordens de magnitude já mencionadas).

ode GertVdE
fonte

Você está executando esta simulação em unidades SI? No mínimo, você deve dividir tudo por algum fator de , para poder eliminar algumas ordens de magnitude.

G * m_{2}

$G*m_{2}$

Oi, eu sou esse, mas ainda não está funcionando ... Os mesmos problemas ocorrem como antes. :(

Você precisa definir sua unidade de massa como uma das massas da lua e suas unidades de comprimento / tempo para definir as coisas como 1. Nada deve ser menor que 1/100 se você a escrever bem. Não há necessidade de um solucionador de balcão. Escreva o código para fazer isso sozinho, onde você controla o tamanho da etapa. Quebras no tamanho da etapa com esses tipos de potenciais podem ocorrer em colisões, onde o solucionador tentará reduzir o tamanho da etapa até a convergência e, na colisão, não há convergência. Você precisa ter certeza de que as órbitas não são colineares; portanto, é necessário visualizar a simulação. Você não pode obter uma resposta como está.

Ron Maimon

Por favor, evite abreviações no título. DGL = Differentialgleichung?

Qual solucionador de ODE padrão você está usando?

Geoff Oxberry 29/03/2012

Respostas:

Sua abordagem atual arruina a estabilidade numérica; na verdade, você provavelmente perde a resolução dessa maneira.

Tome como coordenadas para cada satélite suas variáveis Kepler e o ângulo do plano que contém a posição do satélite, a velocidade e a origem. As equações diferenciais na ausência de interação entre os satélites são então trivialmente simples, e somente a interação se torna um pouco complicada. Como a interação é pequena se os satélites estiverem distantes, a dinâmica resultante deve ser numericamente estável.

Arnold Neumaier
fonte

Em vez de usar um solucionador de ODE "clássico" (rígido), você pode usar algoritmos dedicados para integração numérica geométrica. Veja, por exemplo, este livro e os códigos GNI que você pode encontrar no site da Ernst Hairer .

GertVdE
fonte

Que tal se você tiver três etapas na sua simulação:

atualize a posição de Janus calculando a força de Janus - Saturno.
atualize a posição de Epimeteu calculando a força de Epimeteu - Saturno.
atualize a posição de Janus e Epimetheus calculando a força de Janus - Epimetheus.

Possivelmente usando timesteps mais finos para o # 3.

Não tenho certeza se isso vai ajudar. Suponho que o verdadeiro problema é que a magnitude da força é diferente no caso lua - lua e lua - Saturno, exceto se as luas estiverem próximas?

Alternativamente:

se as luas se fecharem, calcule uma força aproximada de luas - Saturno usando seu vetor de centro de massa e atualize ambas as posições com o mesmo vetor.
se estiverem distantes, atualize-os separadamente.
como antes.

Boa sorte!

fonte