sysROT2={-(Ω^2) X[t]-((d G m M)/((m+M) (((((((d M-(m+M) X[t]))^2)/(((m+M))^2))+((Y[t])^2)))^(3/2))))+((G m X[t])/(((((((d M-(m+M) X[t]))^2)/(((m+M))^2))+((Y[t])^2)))^(3/2)))+((d G m M)/((m+M) (((((((d m+(m+M) X[t]))^2)/(((m+M))^2))+((Y[t])^2)))^(3/2))))+((G M X[t])/(((((((d m+(m+M) X[t]))^2)/(((m+M))^2))+((Y[t])^2)))^(3/2)))-2 Ω Y'[t]+X''[t],-(Ω^2) Y[t]+((G m Y[t])/(((((((d M-(m+M) X[t]))^2)/(((m+M))^2))+((Y[t])^2)))^(3/2)))+((G M Y[t])/(((((((d m+(m+M) X[t]))^2)/(((m+M))^2))+((Y[t])^2)))^(3/2)))+2 Ω X'[t]+Y''[t]}; (* Wyprowadzenie sysROT2 w notatniku 1 *)

W ograniczonym problemie 3 ciał istnieje zachowana ,,energia’’ (całka Jacobiego), którą można zapisać jako sumę energii kinetycznej, odśrodkowej i grawitacyjnej potencjalnej dla trzeciego ciała.

energy=1/2 X'[t]^2+1/2Y'[t]^2-1/2 (X[t]^2+Y[t]^2) Ω^2- G m /Sqrt[(((M/(m+M)d-X[t]))^2)+((Y[t])^2)]-G M/Sqrt[(((m/(m+M)d+X[t]))^2)+((Y[t])^2)]

Mnożąc pierwsze równanie przez X’[t] a drugie przez Y’[t]  i  dodając stronami, wyraz pochodzący od siły Coriolisa kasuje się i równanie można scałkować 1 raz po czasie.

denergy1=sysROT2[[1]]X'[t]+sysROT2[[2]]Y'[t]//Simplify//Expand

Mathematica nie potrafi wykonać takiego całkowania:

(*(Integrate[#,t]&)/@denergy1*)

(*Integrate[denergy1,t]*)

Pozostaje sprawdzić wynik przez różniczkowanie.

denergy2=D[energy,t]

denergy1-denergy2//Expand//Simplify

Rozważmy potencjał U:

U=energy/.X->(x&)/.Y->(y&)

{m=1,M=81,d=1,G=1,Ω=Sqrt[G(m+M)/d^3]} (* Jednostki dla układu Ziemia-Księżyc *)

Widzimy 2 studnie potencjału. Głębsza, pochodzi od ,,Ziemi’’, mniejsza od ,,Księżyca’’. Jest to zgodne z intuicją, że bardzo blisko danego ciała, jego grawitacja musi dominować.

Plot3D[U,{x,-2d,2d},{y,-2d,2d},PlotPoints->64]

Aby powiedzieć coś więcej na temat ruchu w obszarze przejściowym, użyjemy punktów Lagrange’a jako skali.

(* Wzory na punkty Lagrange'a wyprowadzone w notatniku 2 *)

L1=Function[q,{Root[1+q-(q^3)-(q^4)+(2+4 q+5 (q^2)+4 (q^3)+2 (q^4)) #1+(1-(q^4)) ((#1)^2)+(1-2 q-6 (q^2)-2 (q^3)+(q^4)) ((#1)^3)+(2+4 q-4 (q^3)-2 (q^4)) ((#1)^4)+(1+4 q+6 (q^2)+4 (q^3)+(q^4)) ((#1)^5)&,1],0}];
L2=Function[q,{Root[-1-q-(q^3)-(q^4)+(-2-4 q+(q^2)+4 (q^3)+2 (q^4)) #1+(-1-6 q-6 (q^2)-2 (q^3)-(q^4)) ((#1)^2)+(1-2 q-6 (q^2)-2 (q^3)+(q^4)) ((#1)^3)+(2+4 q-4 (q^3)-2 (q^4)) ((#1)^4)+(1+4 q+6 (q^2)+4 (q^3)+(q^4)) ((#1)^5)&,1],0}];
L3=Function[q,{Root[1+q+(q^3)+(q^4)+(2+4 q+(q^2)-4 (q^3)-2 (q^4)) #1+(1+2 q+6 (q^2)+6 (q^3)+(q^4)) ((#1)^2)+(1-2 q-6 (q^2)-2 (q^3)+(q^4)) ((#1)^3)+(2+4 q-4 (q^3)-2 (q^4)) ((#1)^4)+(1+4 q+6 (q^2)+4 (q^3)+(q^4)) ((#1)^5)&,1],0}];
L4 = Function[q,{-(1-q)/(1+q)/2,Sqrt[3]/2}];
L5 = Function[q,{-(1-q)/(1+q)/2,-Sqrt[3]/2}];

Obliczamy wartość potencjału U w punktach Lagrange’a.

energyL1=U/.{x->L1[M/m][[1]],y->0}//N
energyL2=U/.{x->L2[M/m][[1]],y->0}//N
energyL3=U/.{x->L3[M/m][[1]],y->0}//N
energyL4=U/.{x->L4[M/m][[1]],y->L4[M/m][[2]]}//N
energyL5=U/.{x->L5[M/m][[1]],y->L5[M/m][[2]]}//N

energie={energyL1,energyL2,energyL3,energyL4,energyL5};

Jak widać z tabelki, “minimum” energii (lepiej: najmniejsza spośród wartości energii) jest związane z (Subscript[L, 1]). Maksymalna: (Subscript[L, 4,5]).

{energie,Table["U(L"<>ToString[i]<>")",{i,1,5}]}//Transpose//Sort//Reverse//TableForm

lagrangePoints=Graphics3D[
{
{PointSize[Large],Red,Point[Join[L1[M/m],{energyL1}]]},
{PointSize[Large],Red,Point[Join[L2[M/m],{energyL2}]]},
{PointSize[Large],Red,Point[Join[L3[M/m],{energyL3}]]},
{PointSize[Large],Red,Point[Join[L4[M/m],{energyL4}]]},
{PointSize[Large],Red,Point[Join[L5[M/m],{energyL5}]]}
}
];

{minE,maxE}=MinMax[energie]

ΔE=Subtract@@MinMax[{energyL1,energyL2,energyL3,energyL4,energyL5}]//Abs

Wykres poniżej pokazuje, że punkty (Subscript[L, 1,2])są puntkami siodłowymi U(X,Y).

Show[
Plot3D[U,{x,-2d,2d},{y,-2d,2d},PlotRange->{minE-ΔE,maxE+ΔE},PlotPoints->128,ClippingStyle->None],
lagrangePoints
]

ΔEL4L3=energyL4-energyL3

Powiększony jeszcze raz wykres wyjaśnia, że (Subscript[L, 3]) jest punktem siodłowym, natomiast (Subscript[L, 3,4]) to maksimum U(X,Y).

Show[
Plot3D[U,{x,-2d,2d},{y,-2d,2d},PlotRange->{energyL3-ΔEL4L3,maxE+ΔEL4L3},PlotPoints->128,ClippingStyle->None],
lagrangePoints
]

Poniżej pokazuję tzw. krzywe zerowej prędkości.

Manipulate[ContourPlot[U==Energy,{x,-2d,2d},{y,-2d,2d},PlotPoints->50,Contours->20,
Epilog->{
{Red, Point[L1[M/m]],PointSize[Large]},
{Red, Point[L2[M/m]],PointSize[Large]},
{Red, Point[L3[M/m]],PointSize[Large]},
{Red, Point[L4[M/m]],PointSize[Large]},
{Red, Point[L5[M/m]],PointSize[Large]},
{Black, Point[{-m/(m+M),0}],PointSize[Large]},
{Black, Point[{M/(m+M),0}],PointSize[Large]},
Text["(Subscript[L, 1])",L1[M/m],{0,1},BaseStyle->FontSize->12],
Text["(Subscript[L, 2])",L2[M/m],{0,1},BaseStyle->FontSize->12],
Text["(Subscript[L, 3])",L3[M/m],{0,1},BaseStyle->FontSize->12],
Text["(Subscript[L, 4])",L4[M/m],{0,1},BaseStyle->FontSize->12],
Text["(Subscript[L, 5])",L5[M/m],{0,1},BaseStyle->FontSize->12]
}
]

,{Energy,{energyL1,energyL2,energyL3,energyL4,energyL5}}]

Ponieważ energia kinetyczna jest nieujemna, wstawienie zerowych prędkości daje nierówność E>U(X,Y) wyznaczającą dozwolony obszar ruchu. Obszar ten nazywamy Regionem Hilla.

Manipulate[RegionPlot[U<Energy,{x,-2d,2d},{y,-2d,2d},PlotPoints->64,PerformanceGoal->"Quality",
Epilog->{
{Red, Point[L1[M/m]],PointSize[Large]},
{Red, Point[L2[M/m]],PointSize[Large]},
{Red, Point[L3[M/m]],PointSize[Large]},
{Red, Point[L4[M/m]],PointSize[Large]},
{Red, Point[L5[M/m]],PointSize[Large]},
{Black, Point[{-m/(m+M),0}],PointSize[Large]},
{Black, Point[{M/(m+M),0}],PointSize[Large]},
Text["(Subscript[L, 1])",L1[M/m],{0,1},BaseStyle->FontSize->12],
Text["(Subscript[L, 2])",L2[M/m],{0,1},BaseStyle->FontSize->12],
Text["(Subscript[L, 3])",L3[M/m],{0,1},BaseStyle->FontSize->12],
Text["(Subscript[L, 4])",L4[M/m],{0,1},BaseStyle->FontSize->12],
Text["(Subscript[L, 5])",L5[M/m],{0,1},BaseStyle->FontSize->12]
}
],
(*continuous slider*)
{{Energy,energyL1,"Energy"},energyL1-ΔE,energyL4+ΔE,0.01},
(*discrete presets that write into the same variable*)Control[{{Energy,energyL1},{energyL1->"(Subscript[L, 1])",energyL2->"(Subscript[L, 2])",energyL3->"(Subscript[L, 3])",energyL4->"(Subscript[L, 45])"}}]
]

Podsumowanie:
	1. Dla E<U((Subscript[L, 1])) mamy 3 rozłączne obszary dozwolone: 
		a) dookoła ciała o masie M
		b) dookoła ciała o masie m
		c) obszar ,,zewnętrzny’’ rozciągający się do ∞
	2. Dla E=U((Subscript[L, 1])) ZVC (zero-velocity curve) stanowi tzw. krytyczną powierzchnię Roche’a układu 2 ciał. Ma ona ogromne znaczenie dla teorii ewolucji ukłądów podwójnych gwiazd, 			gdyż jest możliwość wypełnienia jednego z ,,lobów’’ i akrecji (transferu masy) przez wewnętrzny punkt (Subscript[L, 1])
	3. Dla U((Subscript[L, 1]))<E<U((Subscript[L, 2])) obszary wokół mas m i M łączą się. Ruch polega na chaotycznym ,,przeskakiwaniu’’ z rejonów m<->M, przedzielonych erami ,,stabilności’’. 
	4. Dla E>U((Subscript[L, 2])) przez zewnętrzny punkt Lagrange’a (Subscript[L, 2]) otwiera się ,,droga ucieczki'' z układu. Wszystkie obszary są połączone. Jest możliwe wyrzucenie ciała z układu do ∞.
	5. Dla E<U((Subscript[L, 45])) istnieją obszary ,,zabronione’ w pobliżu punktów (Subscript[L, 345]). 
	6. Dla E>U((Subscript[L, 45])) ruch jest dozwolony wszędzie, a Region Hilla staje się całą płaszczyzną XY.

Poniżej przykład rozwiązania numerycznego równań ruchu.

(* Warunki początkowe dla połączonych przez L1 regionów mas m i M *)
v0x=-1.3;
v0y=-0.02;
{x0,y0}={0.75M/(m+M)d,-0.2};

(* Warunki początkowe demonstrujace stabilnosc punktu L5 *)
(*
v0x=-0.03;
v0y=-0.02;
{x0,y0}={(M-m)/(m+M)d/2,d Sqrt[3]/2}
*)

U0=U/.{x->x0,y->y0}

E0=1/2  (v0x^2+v0y^2)+U0//N

E0<energyL3

E0<energyL2

E0<energyL1

E0<energyL4

energyL2

tMAX=60

rozw=NDSolve[{sysROT2[[1]]==0,sysROT2[[2]]==0,X[0]==x0,X'[0]==v0x,Y[0]==y0,Y'[0]==v0y},{X,Y},{t,0,tMAX}(*,MaxSteps->10^6,Method->"ExplicitRungeKutta"*)][[1]]

Plot[{energy/.rozw},{t,0,tMAX},GridLines->{None,energie}]

traj=ParametricPlot[{X[τ],Y[τ]}/.rozw,{τ,tMAX*0,tMAX},PlotPoints->128,PlotStyle->Thin]

hillregion=RegionPlot[U<E0,{x,-1.5d,1.5d},{y,-1.5d,1.5d},PlotPoints->128]

w=Show[traj,hillregion,AxesOrigin->{0,0},
Epilog->{Disk[{-m/(m+M)d,0},0.04],Disk[{M/(m+M)d,0},0.01],{Red,Point[L1[m/M]],Point[L2[m/M]],Point[L3[m/M]],Point[-L4[m/M]],Point[-L5[m/M]]}},ImageSize->300]

trajROT=ParametricPlot[({X[τ],Y[τ]}.{{Cos[Ω τ],Sin[Ω τ]},{-Sin[Ω τ],Cos[Ω τ]}})/.rozw,{τ,tMAX/2,tMAX},PlotPoints->128]