Problem treh teles

Literatura :  Vladimir Igorevich Arnold, V. V. Kozlov, A. I. Neĭshtadt:  Mathematical aspects of classical and celestial mechanics.

Začnimo malce bloj splošno in na začetku dovolimo poljubno število teles.

n=3;
MyColourPalete ={Red,Green,Blue,Black};

n je število točk. Za n > 4 je potrebno padaljšati niz MyColourPalete z dodatnimi barvami.

Do[p[k]= {x[k][t],y[k][t]},{k,n}]

Koordinate točk. Omejili se bomo na ravninske primer.

MyNorm[p_]:= (Sqrt[Sum[(((p[[k]]))^2),{k,2}]])

eq[k_]:=Thread[ D[p[k],{t,2}]== Sum[If[i≠ k,((-m[i])/((MyNorm[p[k]-p[i]])^3))(p[k]-p[i]),0],{i,3}]]

Gibalne enačbe v brezdimenzijski obliki.

eqeq=Flatten[Table[eq[k],{k,n}]]

Homotetična rešitev

V začetnem trenutku so telesa v  ogliščih pravilnega mnogokotnika in imajo vsa enako obodno hitrost. Znano je, da se vsako telo giblje po elipsi tako, da vseskozi telesa tvorijo pravilni mnogokotnik. Preverimo to z numeričnim izračunom.

zp =Flatten[ Table[{x[k+1][0]== Cos[(k/(n ))2π],y[k+1][0]== Sin[(k/(n ))2π]},{k,0,n-1}]]
zv= Flatten[  Table[{x[k+1]'[0]== -v0 Sin[(k/(n ))2π],y[k+1]'[0]== v0 Cos[(k/(n ))2π]},{k,0,n-1}]]

tmax=4;
mdata = Join[Table[m[k]-> 1,{k,n}],{v0-> (1/4)}]
eqeqd= Join[eqeq,zp,zv]/.mdata;
slv=NDSolve[eqeqd,Join[Table[x[k],{k,n}],Table[y[k],{k,n}]],{t,0,tmax}];
tbl=Table[{x[k][t],y[k][t]}/.slv[[1]],{k,n}];
slika=ParametricPlot[tbl,{t,0,tmax},PlotStyle-> MyColourPalete];
points=Show[Table[Table[Graphics[{ MyColourPalete[[k]],PointSize[Large],Point[{x[k][t],y[k][t]}/.slv]}],{t,0,tmax,1}],{k,1,n}]];
text     =Show[Table[Table[Graphics[Text[t,{x[k][t],y[k][t]}/.slv[[1]],{0,-1}]],{t,0,tmax,1}],{k,1,n}]];
Show[{slika,points,text}]

Številke nad točkami označujejo čas.

Če podaljšamo računski  interval, dobimo:

tmax=30;
mdata = Join[Table[m[k]-> 1,{k,n}],{v0-> (1/4)}]
eqeqd= Join[eqeq,zp,zv]/.mdata;
slv=NDSolve[eqeqd,Join[Table[x[k],{k,n}],Table[y[k],{k,n}]],{t,0,tmax}];
tbl=Table[{x[k][t],y[k][t]}/.slv[[1]],{k,n}];
slika=ParametricPlot[tbl,{t,0,tmax},PlotStyle-> MyColourPalete];
points=Show[Table[Table[Graphics[{ MyColourPalete[[k]],PointSize[Large],Point[{x[k][t],y[k][t]}/.slv]}],{t,0,tmax,1}],{k,1,n}]];
text     =Show[Table[Table[Graphics[Text[t,{x[k][t],y[k][t]}/.slv[[1]],{0,-1}]],{t,0,tmax,1}],{k,1,n}]];
Show[{slika,points,text}]

Vidimo, da rešitev numerično podivija. Če želimo ohraniti pravo rešitev potrebujemo večjo natančnost računanja.

tmax=30;
mdata = Join[Table[m[k]-> 1,{k,n}],{v0-> (1/4),ag-> 20,pg-> 20,wp->30}]
eqeqd= Join[eqeq,zp,zv]/.mdata;
slv=NDSolve[eqeqd,Join[Table[x[k],{k,n}],Table[y[k],{k,n}]],{t,0,tmax},AccuracyGoal-> ag/.mdata,PrecisionGoal->pg/.mdata,WorkingPrecision->wp/.mdata,MaxSteps->Infinity];
tbl=Table[{x[k][t],y[k][t]}/.slv[[1]],{k,n}];
slika=ParametricPlot[tbl,{t,0,tmax},PlotStyle-> MyColourPalete];
points=Show[Table[Table[Graphics[{ MyColourPalete[[k]],PointSize[Large],Point[{x[k][t],y[k][t]}/.slv]}],{t,0,tmax,1}],{k,1,n}]];
text     =Show[Table[Table[Graphics[Text[t,{x[k][t],y[k][t]}/.slv[[1]],{0,-1}]],{t,0,tmax,1}],{k,1,n}]];
Show[{slika,points,text}]

Rešitev je sedaj prava. Iz slike razberemo tudi brzino gibanja.

Pitagorejska trojica

Poglejmo še vpliv numerične napake pri naslednjem primeru.

zp =Flatten[ {{x[1][0]== 1,y[1][0]==3},{x[2][0]==-2,y[2][0]==-1},{x[3][0]== 1,y[3][0]==-1}}];
zv= Flatten[  Table[{x[k+1]'[0]== 0,y[k+1]'[0]==0},{k,0,n-1}]];
mdata = {m[1]-> 3,m[2]-> 4,m[3]-> 5}

Mase sestavljajo Pitagorejsko trojico, in na začetku mirujejo.

Oglejmo si prvo začetek gibanja

tmax=10;tmin=0;
eqeqd= Join[eqeq,zp,zv]/.mdata;
slv=NDSolve[eqeqd,Join[Table[x[k],{k,n}],Table[y[k],{k,n}]],{t,0,tmax},MaxSteps->20000];
tbl=Table[{x[k][t],y[k][t]}/.slv[[1]],{k,n}];
slika=ParametricPlot[tbl,{t,tmin,tmax},PlotStyle-> MyColourPalete];
points=Show[Table[Table[Graphics[{ MyColourPalete[[k]],PointSize[Large],Point[{x[k][t],y[k][t]}/.slv]}],{t,tmin,tmax,1}],{k,1,n}]];
text     =Show[Table[Table[Graphics[Text[t,{x[k][t],y[k][t]}/.slv[[1]],{0,-1}]],{t,tmin,tmax,1}],{k,1,n}]];
Show[{slika,points,text}]

Sedaj pa za daljši čas.

tmax=70;tmin=0;
eqeqd= Join[eqeq,zp,zv]/.mdata;
slv=NDSolve[eqeqd,Join[Table[x[k],{k,n}],Table[y[k],{k,n}]],{t,0,tmax},MaxSteps->30000];
tbl=Table[{x[k][t],y[k][t]}/.slv[[1]],{k,n}];
slika=ParametricPlot[tbl,{t,tmin,tmax},PlotStyle-> MyColourPalete];
points=Show[Table[Table[Graphics[{ MyColourPalete[[k]],PointSize[Large],Point[{x[k][t],y[k][t]}/.slv]}],{t,tmin,tmax,1}],{k,1,n}]];
text     =Show[Table[Table[Graphics[Text[t,{x[k][t],y[k][t]}/.slv[[1]],{0,-1}]],{t,tmin,tmax,1}],{k,1,n}]];
Show[{slika,points,text}]

Vidimo, da rdeče in zeleno telo v paru ubežita v neskončnost v eno smer, v nasprotno pa dvakrat hitreje modro telo. Ali je rešitev prava? Ponovimo račun z večjo natančnostjo.

tmax=70;
mdata = Join[ {m[1]-> 3,m[2]-> 4,m[3]-> 5},{ag-> 20,pg-> 20,wp->30}]
eqeqd= Join[eqeq,zp ,zv]/.mdata;
slv=NDSolve[eqeqd,Join[Table[x[k],{k,n}],Table[y[k],{k,n}]],{t,0,tmax},AccuracyGoal-> ag/.mdata,PrecisionGoal->pg/.mdata,WorkingPrecision->wp/.mdata,MaxSteps->Infinity];
tbl=Table[{x[k][t],y[k][t]}/.slv[[1]],{k,n}];
slika=ParametricPlot[tbl,{t,0,tmax},PlotStyle-> MyColourPalete];
points=Show[Table[Table[Graphics[{ MyColourPalete[[k]],PointSize[Large],Point[{x[k][t],y[k][t]}/.slv]}],{t,0,tmax,1}],{k,1,n}]];
text     =Show[Table[Table[Graphics[Text[t,{x[k][t],y[k][t]}/.slv[[1]],{0,-1}]],{t,0,tmax,1}],{k,1,n}]];
Show[{slika,points,text}]

Dobimo povsem drugo sliko. Ali je prava?