Simulando um sistema de vários corpos

Um mini sistema solar

O famoso problema dos três corpos (ver por exemplo o link) não possui solução analítica, isto é, não possível sob condições arbitrárias se calcular exatamente o comportamento do sistema para um instante qualquer. Entretanto, é fácil escrever um código (programa) para se obter uma solução numérica aproximada.

Considerando somente o caso plano e as equações de Newton, o código Scilab abaixo pode gerar um mini sistema planetário bem interessante. Dependendo das condições iniciais, é possível ver "planetas" sendo expulsos, "cometas" sendo capturados em uma órbita próxima ao centro ("Sol") como mostra a figura acima, etc.

Código abaixo (quase sem comentários) pode ser expandido para 3 dimensões (eixos x, y e z, mais o tempo) e otimizado para rodar mais rápido.
=======================================
clc;
close;

// sistema (plano) de varios corpos
G = 6.67e-11;
massa = [1500, 2, 1, 1.5, 0.01]*1e6; // massa
px = [0, 1, 2, -3.2, -0.5]; // pos. x
py = [0, 0, 1, 0, 11]; // pos. y

Nmax = max(size(massa));
tmax = 60000;

// velocidades:
vx = 0*px;
vy = [0, 0.25, 0.25, -0.2, -0.1];

vpx = zeros(Nmax,tmax);
vpy = vpx;

dt = 0.005;
Nmax = max(size(massa));

// forças e aceleraçoes:
disp("px, py");
for t=1:tmax
for n=1:Nmax
    ftx = 0;
    fty = ftx;   
    for m=1:Nmax
        if abs(n-m)>0 then
            dx = px(m)-px(n);
            dy = py(m)-py(n);
            d2 = dx*dx + dy*dy;
            fmn = G*massa(m)/d2;
            angnm = atan(dy,dx);
            ftx = ftx + fmn*cos(angnm);
            fty = fty + fmn*sin(angnm);
        end     
    end
        ax = ftx;
        ay = fty;
        vx(n) = vx(n) + ax*dt;
        vy(n) = vy(n) + ay*dt;
        px(n) = px(n) + vx(n)*dt;
        py(n) = py(n) + vy(n)*dt; 
        vpx(n,t)=px(n);
        vpy(n,t)=py(n);
end
end;
disp('pronto');

close;
cores=["b","m","g","k","y"];
for k=Nmax:-1:1
    plot(vpx(k,:), vpy(k,:),cores(k));
end;
legend('C-0.01','P-1.5','P-1','P-2','Sol');
=======================================================