Um pouco de caos no plano: sistemas de três corpos

 

Sistemas com 3 ou mais corpos ligados gravitacionalmente não possuem solução analítica. As leis de Newton são válidas, mas, de forma geral, encontrar uma solução fechada para os movimentos é impossível.  Na verdade, se os três corpos apresentarem massas da mesma ordem de grandeza, a tendência é que o movimento seja caótico e imprevisível a longo prazo.  

Um sistema caótico é um sistema dinâmico que exibe comportamento altamente sensível às condições iniciais, onde pequenas mudanças podem levar a resultados drasticamente diferentes (efeito borboleta). Apesar de serem regidos por regras determinísticas, esses sistemas parecem aleatórios e são imprevisíveis a longo prazo. Exemplos incluem o clima, turbulência de fluidos e certos sistemas biológicos ou econômicos. Características principais incluem não linearidade, dependência sensível às condições iniciais e padrões complexos, como atratores estranhos. Um sistema de três (ou mais) corpos de massas semelhantes interagindo gravitacionalmente é um exemplo clássico de sistema caótico. No problema dos três corpos, as interações gravitacionais mútuas criam um sistema dinâmico não linear, onde pequenas variações nas condições iniciais (como posições ou velocidades) podem levar a trajetórias completamente diferentes ao longo do tempo. Exceto em casos muito específicos (como configurações altamente simétricas), as órbitas não seguem padrões periódicos simples e exibem comportamento caótico, com sensibilidade às condições iniciais e imprevisibilidade a longo prazo. Isso contrasta com o problema de dois corpos, que é geralmente previsível e possui soluções analíticas estáveis. 

Nesta postagem apresentamos um código Scilab que simula um sistema de três corpos com massas da mesma ordem de grandeza. Mudamos as condições iniciais (posições, massas ou velocidades) e observamos resultados completamente distintos. Em alguns casos, um dos corpos pode ser expulso do sistema. 

Código Scilab (feito com auxílio do ChatGPT -5): 

Prompt:  

Vamos simular o problema dos três corpos no plano? Gere um código Scilab que simule o problema dos três corpos no plano, isto é, os três corpos estão se movendo um espaço 2D. Use as equações básicas de gravitação de Newton. Os parâmetros de simulação (massa dos três corpos, posição inicial, velocidades iniciais) devem ficar fáceis de serem editados. A saída será na forma de um gráfico.

...

// =================================================
// Problema dos Três Corpos no plano (2D) - Gravitação de Newton
// Integração RK4 (passo fixo). Saída: gráfico das trajetórias.
// Editar parâmetros NA SEÇÃO "PARÂMETROS DE SIMULAÇÃO".
// =================================================

// ---------- PARÂMETROS DE SIMULAÇÃO (edite aqui) ------------
G      = 1.0;          // Constante gravitacional (use 1 para unidades normalizadas)
m      = [2; 4; 20];    // Massas [m1; m2; m3]
soft   = 0.001;         // "softening" para evitar singularidades (>= 0)
dt     = 0.002;        // Passo de tempo
Tend   = 80;           // Tempo final da simulação
saveEvery = 5;         // Salvar a cada N passos (para acelerar plots/memória)

// Condições iniciais (posições e velocidades em 2D)
// r = [x1 y1; x2 y2; x3 y3]
r0 = [ 2.0  0.0;
       3.0  1.5;
       0.0  0.0];

// v = [vx1 vy1; vx2 vy2; vx3 vy3]
// Sugestão inicial: dar velocidades opostas às duas massas externas
v0 = [ 0.0  5.00;
       0.0  3.00;
       0.0  0.00];

// Opcional: forçar o centro de massa em repouso (zera momento total)
Pcm = m(1)*v0(1,:) + m(2)*v0(2,:) + m(3)*v0(3,:);
Mtot = sum(m);
v0 = [v0(1,:) - Pcm/Mtot;
      v0(2,:) - Pcm/Mtot;
      v0(3,:) - Pcm/Mtot];

// ------------- FUNÇÕES AUXILIARES ---------------------------
function a=aceleracoes(r, m, G, soft)
    // r: 3x2 (linhas = corpos, colunas = x,y)
    // a: 3x2
    a = zeros(3,2);
    for i = 1:3
        ai = [0, 0];
        for j = 1:3
            if j <> i then
                dx = r(j,1) - r(i,1);
                dy = r(j,2) - r(i,2);
                d2 = dx*dx + dy*dy + soft*soft; // softening
                inv3 = 1 / (d2*sqrt(d2));        // 1 / |r|^3
                ai = ai + G * m(j) * [dx, dy] * inv3;
            end
        end
        a(i,:) = ai;
    end
endfunction

function [rnext, vnext]=rk4_step(r, v, m, G, soft, dt)
    // k1
    a1 = aceleracoes(r, m, G, soft);
    k1r = v;
    k1v = a1;

    // k2
    r2 = r + 0.5*dt*k1r;
    v2 = v + 0.5*dt*k1v;
    a2 = aceleracoes(r2, m, G, soft);
    k2r = v2;
    k2v = a2;

    // k3
    r3 = r + 0.5*dt*k2r;
    v3 = v + 0.5*dt*k2v;
    a3 = aceleracoes(r3, m, G, soft);
    k3r = v3;
    k3v = a3;

    // k4
    r4 = r + dt*k3r;
    v4 = v + dt*k3v;
    a4 = aceleracoes(r4, m, G, soft);
    k4r = v4;
    k4v = a4;

    rnext = r + dt*(k1r + 2*k2r + 2*k3r + k4r)/6;
    vnext = v + dt*(k1v + 2*k2v + 2*k3v + k4v)/6;
endfunction

function [K, U, E]=energia_total(r, v, m, G)
    // Cinética
    K = 0;
    for i=1:3
        vi2 = v(i,1)^2 + v(i,2)^2;
        K = K + 0.5*m(i)*vi2;
    end
    // Potencial gravitacional
    U = 0;
    for i=1:3
        for j=i+1:3
            dx = r(j,1) - r(i,1);
            dy = r(j,2) - r(i,2);
            rij = sqrt(dx*dx + dy*dy);
            U = U - G*m(i)*m(j) / rij;
        end
    end
    E = K + U;
endfunction

clc; /// limpar o 'console'
close(winsid()); /// fecha todas as janelas

// ------------- MALHA DE TEMPO E ARRAYS ----------------------
Nsteps = round(Tend/dt);
Nsaved = floor(Nsteps/saveEvery) + 1;

R1 = zeros(Nsaved, 2); R2 = zeros(Nsaved, 2); R3 = zeros(Nsaved, 2);
times = zeros(Nsaved, 1);
Ehist = zeros(Nsaved, 1);

// Estados
r = r0;  v = v0;

k = 1;
R1(k,:) = r(1,:); R2(k,:) = r(2,:); R3(k,:) = r(3,:);
times(k) = 0;
//[~,~,Ehist(k)] 
r_et = energia_total(r, v, m, G);

// ---------------- INTEGRAÇÃO PRINCIPAL ----------------------
for n = 1:Nsteps
    [r, v] = rk4_step(r, v, m, G, soft, dt);

    if modulo(n, saveEvery) == 0 then
        k = k + 1;
        R1(k,:) = r(1,:); R2(k,:) = r(2,:); R3(k,:) = r(3,:);
        times(k) = n*dt;
        // [~,~,Ehist(k)] 
        r_et = energia_total(r, v, m, G);
    end
end

// ------------------- GRÁFICOS -------------------------------
scf(1); clf(); // Trajetórias
plot(R1(:,1), R1(:,2), "r-");
plot(R2(:,1), R2(:,2), "m-");
plot(R3(:,1), R3(:,2), "k-");
fim = max(size(R1));
// posições finais:
plot([R1(fim,1) R2(fim,1) R3(fim,1)],[R1(fim,2) R2(fim,2) R3(fim,2)],"o"); 
// posições iniciais:
plot([r0(1,1) r0(2,1) r0(3,1)], [r0(1,2) r0(2,2) r0(3,2)], "x"); 
xgrid();
// axis("equal");
xlabel("x"); ylabel("y");
legend(["Corpo 1"; "Corpo 2"; "Corpo 3"; "Posição final"; "Posição inicial"],-1);
title(msprintf("Três Corpos 2D (G=%.3f, dt=%.4f, soft=%.3f)", G, dt, soft));

Exemplos de resultados: