Poniżej wyprowadzimy 5 statycznych rozwiązań RCPTBP (Restricted Circular Planar Three Body Problem) znanych szerzej jako punkty Lagrange’a (Subscript[L, 1])-(Subscript[L, 5]). Ich istnienie jest znaną, ale nieco zaskakującą cechą problemu 3 ciał.

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]==0,-(Ω^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]==0}; (* Wyprowadzenie ukł. równ. w notatniku 1 *)

eqX=First[sysROT2[[1]]]

eqY=First[sysROT2[[2]]]

Istnienie aż 5 niezależnych od czasu rozwiązań powyższego układu równań wydaje się mało oczywiste. Trzy z nich ((Subscript[L, 1])-(Subscript[L, 3])) spodziewamy się znaleźć na linii łączącej ciała m i M (Y=0), jako że GDZIEŚ siły grawitacji i odśrodkowa powinny się równoważyć. Pozostałe ((Subscript[L, 4, ])(Subscript[L, 5])) układają się na wierzchołkach trójkąta równobocznego (!).

wTriangle=Block[{M=3,m=1,G=1,d=1,Ω=Sqrt[G (m+M)/d^3],T=2 Pi/Ω,
(*Pre-calculate positions for clarity*)
posM=M/(m+M) {Cos[Ω T],Sin[Ω T]},
posm=-m/(m+M) {Cos[Ω T],Sin[Ω T]},
posL4={d (M-m)/(M+m)/2,d Sqrt[3]/2},
posL5={d (M-m)/(M+m)/2,-d Sqrt[3]/2}},
ParametricPlot[{0,0},{t,0,T},
Epilog->{Disk[posM,0.02],Disk[posm,0.05]},PlotStyle->{{Red,Thick},{Blue,Thick}},AxesLabel->{X,Y},PlotRange->All,Prolog->{(*Draw Polygons*)
{EdgeForm[{Black,Thin}],FaceForm[White],Polygon[{posM,posm,posL4}]},
{EdgeForm[{Black,Thin}],FaceForm[White],Polygon[{posM,posm,posL5}]},
(*Draw Points*){Red,PointSize[Large],Point[posL5]},
{Blue,PointSize[Large],Point[posL4]},
(*Add Text Labels*)
{Black,
Text[Style["M",Medium],posM,{0,-1}],
Text[Style["m",Medium],posm,{0,-2}],
Text[Style["L4",Medium],posL4,{-2,0}],
Text[Style["L5",Medium],posL5,{2,0}]}}]
]

{eqX,eqY}/.{X->(d(M-m)/(M+m)/2&), Y->( -d Sqrt[3]/2&),Ω->Sqrt[G(m+M)/d^3]} (* współrzędne punktów L4 i L5, trójkąt równoboczny *)

Simplify[%,d>0]

Znalezienie punktów kolinearnych jest kłopotliwe. Można szukać warunków równoważenia się sił grawitacyjnych i odśrodkowej. My spróbujemy dojść do tego wprost z równań.

{eqX,eqY}/.{Y->(0&),X->(X0&)}

rownanie1=%[[1]]/. Ω->Sqrt[G (M+m)/d^3]

rownanie2 = rownanie1 /.{m->q M}//Simplify (* Wprowadzamy stosunek mass q, mass ratio *)

rownanie3=Collect[rownanie2,M]

rownanie4=Last[rownanie3]

rownanie5 = Simplify[rownanie4/.X0->x d,d>0] (* Używamy odległości pomiędzy masami m i M, czyli d, jako jednostki sprowadzającej problem do postaci bezwymiarowej *)

rownanie6=FullSimplify[Numerator[rownanie5],q>0&&x∈Reals&&d>0]

rownanie7 = Last[rownanie6](* Czemu Mathematica NIE WIE że z/Abs[z]^3 = 1/z/Abs[z] ? ????    *)

Simplify[z/Abs[z]^3==1/z/Abs[z],z∈Reals]

FullSimplify[rownanie7, x∈Reals&&q>0]

Solve[rownanie7==0&&q>0,x,Reals]

Nieco ładniejsza ręcznie uproszczona wersja równania.

Solve[x==(1+q) ((q/((1+x+x q)Abs[(1+x+x q)]))+(1/((x-q+x q)Abs[(x-q+x q)])))&&q>0,x,Reals]

Refine[%,q>0]

Zapisujemy pozycje kolinearnych punktów (Subscript[L, 1]), (Subscript[L, 2]) i (Subscript[L, 3]).

L123=x/.%;

Położenie kolinearnych punktów Lagrange’a oraz mas m,M w zależności od stosunku mas q=m/M.

Plot[{L123,-1,1,q/(1+q),-1/(1+q)},{q,0,2},PlotRange->{-2,2},PlotStyle->{{Red,Thick},{Green,Thick},{Blue,Thick},{Black,DotDashed,Thin},{Black,DotDashed,Thin},{Black,Thick},{Black,Thin},{Black,Thick}},AspectRatio->1,AxesLabel->{"m/M","r/R"},Epilog->{
Text["(Subscript[L, 1])",{2,L123[[3]]/.q->2},{1,1},BaseStyle->{FontSize->22}],
Text["(Subscript[L, 2])",{2,L123[[1]]/.q->2},{1,-1},BaseStyle->{FontSize->22}],
Text["(Subscript[L, 3])",{2,L123[[2]]/.q->2},{1,-1},BaseStyle->{FontSize->22}],
Text["m",{1,-1/2},{1,-1},BaseStyle->{FontSize->16}],
Text["M",{1,1/2},{1,-1},BaseStyle->{FontSize->16}]
},BaseStyle->{FontSize->16}]

Wyniki powinne być symetryczne ze względu na zamianę m <-> M. Widać to ładnie na wykresie półlogarytmicznym.

LogLinearPlot[{L123,-1,1,q/(1+q),-1/(1+q)},{q,10^-5,10^5},PlotRange->{-3/2,3/2},PlotStyle->{{Red,Thick},{Green,Thick},{Blue,Thick},{Black,DotDashed,Thin},{Black,DotDashed,Thin},{Black,Thick},{Black,Thin},{Black,Thick}},AspectRatio->1,AxesLabel->{"m/M","r/R"},Epilog->{
Text["(Subscript[L, 2])",{2,L123[[3]]/.q->2},{1,1},BaseStyle->{FontSize->22}],
Text["(Subscript[L, 1])",{2,L123[[1]]/.q->2},{1,-1},BaseStyle->{FontSize->22}],
Text["(Subscript[L, 3])",{2,L123[[2]]/.q->2},{1,-1},BaseStyle->{FontSize->22}],
Text["m",{1,-1/2},{1,-1},BaseStyle->{FontSize->16}],
Text["M",{1,1/2},{1,-1},BaseStyle->{FontSize->16}]
},BaseStyle->{FontSize->16}]

Zwizualizujemy punkty Lagrange’a.

L123[[1]]

L1=Function[q,{Evaluate[L123[[1]]],0}//Evaluate]

L2=Function[q,{Evaluate[L123[[2]]],0}//Evaluate]

L3=Function[q,{Evaluate[L123[[3]]],0}//Evaluate]

L4 = Function[q,{-(1-q)/(1+q)/2,Sqrt[3]/2}];

L5 = Function[q,{-(1-q)/(1+q)/2,-Sqrt[3]/2}];

m=Function[q,{-1/(1+q),0}];

M=Function[q,{q/(1+q),0}];

Mtot=0.1

Manipulate[Graphics[
{
Disk[M[q], Mtot/(1+q)^(1/3)],
Disk[m[q],(q/(1+q))^(1/3) Mtot],
{Green,PointSize[0.02],
Point[L1[q]],
Point[L2[q]],
Point[L3[q]],
Point[L4[q]],
Point[L5[q]]},
{Dashed,Blue,Line[{M[q],m[q],L4[q],M[q],L5[q],m[q],L3[q],L1[q]}]},
{
Text["(Subscript[L, 1])",L1[q],{1,1},BaseStyle->FontSize->24],
Text["(Subscript[L, 2])",L2[q],{1,1},BaseStyle->FontSize->24],
Text["(Subscript[L, 3])",L3[q],{-1,1},BaseStyle->FontSize->24],
Text["(Subscript[L, 4])",L4[q],{1,1},BaseStyle->FontSize->24],
Text["(Subscript[L, 5])",L5[q],{1,1},BaseStyle->FontSize->24],
Text["m",M[q],{0,-1},BaseStyle->FontSize->24],
Text["M",m[q],{0,-2},BaseStyle->FontSize->24]
},
Arrow[3/2{{0,-1},{0,1}}],
Arrow[3/2{{-1,0},{1,0}}]
},PlotRange->3/2{{-1,1},{-1,1}}],{{q,4},1,232}]