Aloittelijan kysymys gravitaatiosta

Jos kappaleita on useampi kuin kaksi, alkaa olla jo vaikeuksia. Täytyy ottaa huomioon, että muut kappaleet vaikuttavat paitsi kappaleeseen X, niin myös toisiinsa. Lisäksi kappale X vaikuttaa muihin kappaleisiin, joka taas muuttaa X:ään kohdistuvaa voimaa. Joudutaan hyvin pian liian monimutkaisten yhtälöiden kanssa tekemisiin.

ei taida 3D riittää

tällaista:
http://en.wikipedia.org/wiki/N-body_problem
Miksiköhän m_j ei ole supistettu puolittain ?

Joten, jos näin on, niin painu vittuun!

Analyyttisesta ratkaisusta on turha haaveillakkaan, mutta tehtävä on melko helppo ratkaista numeerisesti.

Tarkoitan siis nyt yksinkertaista pistemäisten kappaleiden ongelmaa kun kappaleiden välillä vaikuttaa gravitaatio. Kyseessähän on vain Newtonin liikeyhtälöiden

F_xi = m_i*x_i´´
F_yi = m_i*y_i´´
F_zi = m_i*z_i´´

ratkaisu jokaiselle kappaleelle i. Voimatermit ovat aukikirjoitettuna kuten aiemmin mainitussa Wikipediassa http://en.wikipedia.org/wiki/N-body_problem .

Kun yhtälöryhmä muutetaan ensimmäisen kertaluvun normaaliryhmäksi, on tehtävässä ratkaistavana 6*n differentiaaliyhtälöä, jossa n on kappaleiden lukumäärä. Kyseessä siis alkuarvotehtävä, jossa kappaleiden alkupaikat ja -nopeudet tunnetaan.

Olen tämmösen joskus aiemmin pyöräyttänyt matlabiin:
----------------------------------------------
function kolme
t = 3;
G = 1;
m = [10; 2; 1];
R0 = [
0 0 0
1 0 0
0 1 0
];
V0 = [
0 0 0
0 1 1
-1 0 1
];

[T,B] = nbody(t,R0,V0,m,G);

function [T,B] = nbody(t,R0,V0,m,G)
n = length(m);
b0 = zeros(6*n,1);
b0(1:(3*n)) = R0(:);
b0((3*n+1):end) = V0(:);
[T,B] = ode15s(@odefun,[0 t],b0,[],m,G);
X = B(:,1:n);
Y = B(:,(n+1):(2*n));
Z = B(:,(2*n+1):(3*n));

plot3(X,Y,Z);
grid on
axis equal
xlabel('x')
ylabel('y')
zlabel('z','rot',0)
legend('kpl 1','kpl 2','kpl 3')

function dbdt = odefun(t,b,m,G)
n = length(m);
k = length(b);
x = b(1:n);
y = b((n+1):(2*n));
z = b((2*n+1):(3*n));
dbdt = zeros(k,1);
apu = zeros(n,3);
for i = 1:(3*n)
dbdt(i) = b(3*n+i);
end
for j = 1:n
apu(j,:) = acc(j,x,y,z,m,G)';
end
dbdt((3*n+1):k) = apu(:);

% Kappaleen i kiihtyvyysvektori
function a_i = acc(i,x,y,z,m,G)
n = length(m);
a_i = zeros(3,1);
j = [1:(i-1),(i+1):n];
apu = G * m(j) ./ ((x(j) - x(i)).^2 + (y(j) - y(i)).^2 + (z(j) - z(i)).^2).^(3/2);
a_i(1) = sum(apu .* (x(j) - x(i)));
a_i(2) = sum(apu .* (y(j) - y(i)));
a_i(3) = sum(apu .* (z(j) - z(i)));
-------------------------------------------------

Eli se ratkaisee kolmen kappaleen ongelman, kun kappaleet ovat aluksi pisteissä
r1 = [0,0,0]
r2 = [1,0,0]
r3 = [0,1,0]

ja alkunopeudet

v1 = [0,0,0]
v2 = [0,1,1]
v3 = [-1,0,1]

Massat ovat

m1 = 10
m2 = 2
m3 = 1

ja yleinen gravitaatiovakio

G = 1

Nuo arvot olivat ihan sattumanvaraisesti valittuja. En rajannut algoritmia vain kolmeen kappaleeseen, vaan siinä voisi olla vaikka 13 kappaletta, mutta selkeyden vuoksi kolme on parempi kuvassa. Tämän näköistä käppyrää tuli:
http://img3.freeimagehosting.net/image.php?ec0aff16a9.jpg

, eli melko kaoottista. Viisaammilla lähtöarvoilla saisi jotain kauniimpaakin.