Visitando novamente o problema dos três corpos.

O Problema dos Três Corpos: quando a física vira caos 

O problema dos três corpos descreve o movimento de três massas pontuais sob a influência mútua da Gravitação Universal. Diferente do problema de dois corpos, não existe uma solução analítica geral em forma fechada para este problema, exigindo métodos numéricos para a integração das trajetórias. Já visitamos esse problema aqui, aqui e aqui.

Na física, poucas coisas são tão simples e tão impossíveis quanto o problema dos três corpos. A pergunta é direta: imagine três objetos no espaço — três estrelas, ou Sol, Terra e Lua -- se atraindo pela gravidade. Dá pra prever exatamente onde cada um vai estar daqui a mil anos? Com dois corpos é fácil. Newton resolveu: órbitas elípticas, tudo certinho e previsível (calculável). Com três, o universo joga as regras no fundo de um abismo (você teria coragem de olhar para o fundo de um abismo?).

Por que ele gera caos?

1. Efeito borboleta gravitacional 
   Mude a posição de um dos corpos em 1 milímetro e, depois de um tempo, todo o sistema pode ter um destino diferente. Um corpo pode ser ejetado pra fora, os outros dois podem colidir, ou tudo pode virar um balé complexo que nunca se repete.

2. Não existe “fórmula mágica” 
   Henri Poincaré provou em 1890: não dá pra escrever uma equação que resolva o problema de uma vez pra sempre. A única saída é simular no computador, passo a passo. E mesmo assim, erros minúsculos no início explodem com o tempo.

3. Ordem no meio do caos  
   Nem todo trio é caótico. Se um corpo for bem levinho, tipo uma nave espacial entre Terra e Lua, o sistema fica previsível. Também existem soluções raras e bonitas, como a órbita em “figura de oito”: três corpos iguais se perseguindo num 8 perfeito no espaço.

O legado.  
O problema dos três corpos não foi “resolvido” -- ele nos ensinou que nem tudo pode ser resolvido. Foi estudando ele que Poincaré criou a teoria do caos, mostrando que sistemas determinísticos podem ser imprevisíveis na prática.

É por isso que sondas espaciais precisam corrigir a rota o tempo todo, e por que não dá pra calcular o destino exato de três estrelas dançando juntas daqui a um milhão de anos.

A física, às vezes, inclui o improviso no roteiro. 

Formulação Matemática

De acordo com a Segunda Lei de Newton e a Lei da Gravitação Universal, a aceleração de cada corpo $i$ é dada pela soma das forças exercidas pelos outros dois corpos $j$ e $k$:

$$\begin{aligned} \frac{d^2 \mathbf{r}_i}{dt^2} &= \sum_{j \neq i} \frac{G m_j (\mathbf{r}_j - \mathbf{r}_i)}{|\mathbf{r}_j - \mathbf{r}_i|^3} \end{aligned}$$ 

Para implementar isso no Scilab, transformamos o sistema de equações diferenciais de segunda ordem em um sistema de primeira ordem. Definimos o vetor de estado $Y$ como:
$$Y = [x_1, y_1, v_{x1}, v_{y1}, x_2, y_2, v_{x2}, v_{y2}, x_3, y_3, v_{x3}, v_{y3}]^T$$ 

Implementação em Scilab

O script abaixo utiliza a função ode para resolver as equações de movimento. Foram escolhidas condições iniciais que geram órbitas complexas e visualmente interessantes.

 clc; close; close;

// Problema dos Tres Corpos em Scilab
// Definicao das constantes e parametros
G = 1; // Constante gravitacional normalizada
m1 = 1.0; m2 = 1.5; m3 = 1.0; // Massas dos corpos

// Funcao que define as derivadas do sistema 
// (Equacoes de Movimento)
function dy=tres_corpos(t, y)
    // Extração de posições e velocidades 
    // do vetor de estado y
    r1 = [y(1); y(2)];
    v1 = [y(3); y(4)];
    r2 = [y(5); y(6)];
    v2 = [y(7); y(8)];
    r3 = [y(9); y(10)];
    v3 = [y(11); y(12)];
    
    // Calculo das distancias relativas ao cubo
    d12 = norm(r2 - r1)^3;
    d13 = norm(r3 - r1)^3;
    d23 = norm(r3 - r2)^3;
    
    // Aceleracoes baseadas na Lei da Gravitacao
    a1 = G * m2 * (r2 - r1) / d12 + G * m3 * (r3 - r1) / d13;
    a2 = G * m1 * (r1 - r2) / d12 + G * m3 * (r3 - r2) / d23;
    a3 = G * m1 * (r1 - r3) / d13 + G * m2 * (r2 - r3) / d23;
    
    // Retorno das derivadas [dr/dt; dv/dt]
    dy = [v1(1); v1(2); a1(1); a1(2); ..
          v2(1); v2(2); a2(1); a2(2); ..
          v3(1); v3(2); a3(1); a3(2)];
endfunction

// Condicoes Iniciais: [x, y, vx, vy] para cada corpo
// Configuracao de Lagrange/Euler modificada para instabilidade
y0 = [ 0.0;  0.0;  0.0;  -0.4;  // Corpo 1
      -0.5;  0.8; -0.5; -0.4;  // Corpo 2
      -0.5; -0.8;  0.5; -0.4]; // Corpo 3

// Definicao do tempo de simulacao
// t = linspace(0, 20, 2000);
t = linspace(0, 20, 2000);

// Resolucao numerica usando o solver 'ode'
sol = ode(y0, 0, t, tres_corpos);

// Extracao das trajetorias para plotagem
plot(sol(1,:), sol(2,:), 'r'); // Trajetoria Corpo 1
plot(sol(5,:), sol(6,:), 'g'); // Trajetoria Corpo 2
plot(sol(9,:), sol(10,:), 'b'); // Trajetoria Corpo 3

// Formatacao do grafico
xtitle("Problema dos Tres Corpos - Trajetorias", ...
"Posicao X", "Posicao Y");
legend(["Corpo 1"; "Corpo 2"; "Corpo 3"]);
xgrid();