Qual é a maneira correta de integrar em simulações de astronomia?

Estou criando um simulador simples de astronomia que deve usar a física newtoniana para simular o movimento de planetas em um sistema (ou qualquer objeto, por sinal). Todos os corpos são círculos em um plano euclidiano, que possuem propriedades como posição, velocidade, massa, raio e a força resultante.

Quero atualizar o universo em pequenos intervalos de tempo, geralmente alguns milissegundos, mas não sei como calcular corretamente as alterações de posição.

A força é simples: fr = sum(G * body.m * bodyi.m / dist(body, bodyi)^2).

Mas como eu continuo a partir daí?

Eu poderia fazer isso:

a = Fr/body.m
v += a*dt
position += v*dt

Mas isso seria, obviamente, falso. Talvez se eu adicionasse 0,5 como um fator no cálculo da posição?

Você basicamente tem a resposta - não é necessário o fator 0,5.

Essencialmente, você possui um sistema bidimensional de EDOs de primeira ordem: onde tudo é função do tempo, exceto presumivelmentem, e pontos indicam derivadas do tempo. Se você fizer uma diferenciação simples de primeira ordem no estilo Euler avançado, encontrará x n + 1 -x

$$\begin{align} \dot{x} & = v \\ \dot{v} & = \frac{F}{m}, \end{align}$$ $m$ ou xn+ $$\begin{align} \frac{x_{n+1}-x_n}{\Delta t} & = v_n \\ \frac{v_{n+1}-v_n}{\Delta t} & = \frac{F_n}{m}, \end{align}$$ $$\begin{align} x_{n+1} & = x_n + \Delta t \cdot v_n \\ v_{n+1} & = v_n + \Delta t \frac{F_n}{m}. \end{align}$$ $n$

$t_n$$t_{n+1}$$t_{n+1/2}$$x_0$$v_{1/2}$

$$\begin{align} x_{n+1} & = x_n + \Delta t \cdot v_{n+1/2} \\ v_{n+1/2} & = v_{n-1/2} + \Delta t \frac{F_n}{m} \end{align}$$ para integrar adiante no tempo. Isso é conhecido como método saltador . Com isso, seu sistema conserva uma espécie de energia e é menos provável que as órbitas voem para o infinito ou algo parecido devido ao crescimento exponencial do erro de arredondamento.

O único problema é como obter $v_{1/2}$, já que provavelmente você começa com $v_0$. Lá, você deve usar um esquema tão preciso quanto você está escrevendo, sendo o Runge-Kutta de quarta ordem uma escolha popular. Pode ser instável a longo prazo, mas há tantos erros que você introduzirá em meio intervalo de tempo, e esse erro será mantido pequeno posteriormente pelo esquema de salto.

Finalmente, essa resposta se aplica a qualquer sim geral de gravidade newtoniana. Se você realmente deseja círculos perfeitos , como mencionado na passagem da pergunta, não os conseguirá, exceto em um sistema idealizado no qual os planetas não interagem entre si e as condições iniciais são escolhidas da maneira certa. Se for esse o caso, você não precisa se integrar, pois a velocidade angular (radianos por unidade de tempo) de um objeto é simplesmente

$$\omega = \sqrt{\frac{GM}{r^3}},$$ Onde $M$ é a massa do objeto central e $r$é o raio da órbita. Isso pode ser usado para testar a precisão da sua simulação.