Pitagorejski problem treh teles

Tri masne točke, z masami 3, 4 in 5 na začetku mirujejo v točkah s koordinatami (1,3), (-2,-1) in (1,-1).
Sile med točkami so gravitacijske sile.
Določi gibanje

vecr[i_,t_] := {x[i,t],y[i,t]}

Krjevni vektorji do točk.

MyNorm[vec_] := (Sqrt[vec.vec])

Gibalne enačbe.

eq1=Thread[m1 D[vecr[1,t],{t,2}]==- ((κ m1 m2 (vecr[1,t]-vecr[2,t]))/((MyNorm[vecr[1,t]-vecr[2,t]])^3))-  ((κ m1 m3 (vecr[1,t]-vecr[3,t]))/((MyNorm[vecr[1,t]-vecr[3,t]])^3))]
eq2=Thread[m2 D[vecr[2,t],{t,2}]==- ((κ m2 m1 (vecr[2,t]-vecr[1,t]))/((MyNorm[vecr[2,t]-vecr[1,t]])^3))-  ((κ m2 m3 (vecr[2,t]-vecr[3,t]))/((MyNorm[vecr[2,t]-vecr[3,t]])^3))]
eq3=Thread[m3 D[vecr[3,t],{t,2}]==- ((κ m3 m1 (vecr[3,t]-vecr[1,t]))/((MyNorm[vecr[3,t]-vecr[1,t]])^3))-  ((κ m3 m2 (vecr[3,t]-vecr[2,t]))/((MyNorm[vecr[3,t]-vecr[2,t]])^3))]

Podatki

podatki = {m1-> 3,m2-> 4,m3-> 5,κ-> 1}

Začetni pogoji

ZacetniPogoji =Flatten[Map[Thread,{vecr[1,0]== {1,3},vecr[2,0]== {-2,-1},vecr[3,0]== {1,-1},D[vecr[1,t],t]== 0,D[vecr[2,t],t]== 0,D[vecr[3,t],t]== 0}]/.t-> 0]

eqeq = Join[Flatten[Join[eq1,eq2,eq3]]/.podatki,ZacetniPogoji]

Sistem diferencialnih enačb rešimo numerično.
Izkaže se, da z numeričnim izračunom ne moremo predvideti, kako se sistem obnaša za velik t.

Izračun pri standardni natančnosti

Tiri točk so obarvani rdeče, zeleno in modro.

tmax =50;
resitev=NDSolve[eqeq,Flatten[Table[vecr[i,t],{i,3}]],{t,0,tmax}];
slika1=ParametricPlot[{vecr[1,t]/.resitev[[1]],vecr[2,t]/.resitev[[1]],vecr[3,t]/.resitev[[1]]},{t,0,tmax},PlotStyle->{Red,Green,Blue}]

tmax =120;
resitev=NDSolve[eqeq,Flatten[Table[vecr[i,t],{i,3}]],{t,0,tmax}];
slika1=ParametricPlot[{vecr[1,t]/.resitev[[1]],vecr[2,t]/.resitev[[1]],vecr[3,t]/.resitev[[1]]},{t,0,tmax},PlotStyle->{Red,Green,Blue}]

V neskončnost uide par rdeči/zeleni.(Rezultat odvisen od verzije Mathematice!)

Izračun pri večji natančnosti

Račun je časovno zahtevnejši.

tmax =160;
resitev=NDSolve[eqeq,Flatten[Table[vecr[i,t],{i,3}]],{t,0,tmax},MaxSteps->50000, WorkingPrecision->22];
slika1=ParametricPlot[{vecr[1,t]/.resitev[[1]],vecr[2,t]/.resitev[[1]],vecr[3,t]/.resitev[[1]]},{t,0,tmax},PlotStyle->{Red,Green,Blue}]

Sedaj uide par zeleno modri.