Poniżej przedstawiam wyprowadzenie równań płaskiego, kołowego, ograniczonego problemu 3 ciał poprzez jawne przejście od układu inercjalnego do korotującego. Dwie masy: m,M rotują z prędkością kątową (Ω^2)=G(m+M)/(d^3) .

Równania ruchu w układzie inercjalnym. Masy wytwarzające pole grawitacyjne poruszają się po okręgach.

eqX=x''[t]+  G M/Sqrt[(x[t]+m/(m+M)d Cos[Ω*t])^2+(y[t]+m/(m+M)d Sin[Ω*t])^2]^3 (x[t]+m/(m+M) d*Cos[Ω*t])+G m (x[t]-M/(m+M)d Cos[Ω*t])/Sqrt[(x[t]-M/(m+M)d Cos[Ω*t])^2+(y[t]-M/(m+M)d Sin[Ω*t])^2]^3

eqY=y''[t]+  G M/Sqrt[(x[t]+m/(m+M)d Cos[Ω*t])^2+(y[t]+m/(m+M)d Sin[Ω*t])^2]^3 (y[t]+m/(m+M) d*Sin[Ω*t])+G m (y[t]-M/(m+M)d Sin[Ω*t])/Sqrt[(x[t]-M/(m+M)d Cos[Ω*t])^2+(y[t]-M/(m+M)d Sin[Ω*t])^2]^3

o={{Cos[Ω t],-Sin[Ω t]},{Sin[Ω t],Cos[Ω t]}};MatrixForm[o] (* Macierz transforamacji do układu rotującego z keplerowską prędkością kątową Ω *)

Transformujemy układ współrzędnych do układu korotującego X,Y.

sysROT={eqX,eqY}/.{x->(X[#]Cos[Ω #]-Y[#]Sin[Ω #]&),y->(X[#]Sin[Ω #]+Y[#]Cos[Ω #]&)}

Dodajemy pomnożone przez sin i cos równania stronami i upraszczamy.

FullSimplify[-sysROT[[1]]Sin[Ω t]+sysROT[[2]]Cos[Ω t],d>0]

Otrzymujemy proste równanie na Y. Widać człony Coriolisa i odśrodkowy.

eqY=Simplify[%,m>0&&M>0&&G>0&&d>0]

Nieco innym (w celach dydaktycznych) sposobem wyprowadzamy równanie na X.

Simplify[sysROT[[1]]/Sin[Ω t]+sysROT[[2]]/Cos[Ω t],d>0]

Expand[%]

Collect[%,Cot[t Ω]]

Collect[%,Tan[Ω t]]

Simplify[%]

%/(2 Csc[2 Ω t])

eqX=FullSimplify[%,M>0&&m>0&&R>0&&G>0]

Otrzymujemy układ równań różniczkowych zwyczajnych II rzędu na niewiadome funkcje X,Y.

sysROT2={eqX==0,eqY==0}

Powyższy układ równań płaskiego (bo ruch odbywa się wyłącznie w płaszczyźnie orbitalnej mas m i M) , kołowego (masy m i M poruszają się po orbitach kołowych z częstością Ω), ograniczonego ( bo masa trzeciego ciała μ->0 ) problemu trzech ciał można napisać od razu, identyfikując siły nieinercjalne:

-Ω^2 r (siła odśrodkowa)
2 Ω x v (siła Coriolisa *)

Wyprowadzenie w układzie korotującym (nieinercjalnym) XY.

R={X[t],Y[t]}; (* Pozycja 3 ciała o masie μ->0 w układzie korotującym XY *)

Rm=d{ M/(m+M),0}; (* Pozycja ciała o masie m, nieruchoma w układzie korotującym *)

RM=d{-m/(m+M),0};(* Pozycja ciała o masie M, nieruchoma w układzie korotującym *)

sysROT3=D[R,{t,2}]-
Ω^2 R (* Siła odśrodkowa *) +
2 Ω Cross[D[R,t]] (* Siła Coriolisa *)+
G m (R-Rm)/((R-Rm).(R-Rm))^(3/2) (* Siła grawitacji newtonowskiej pochodząca od ciała o masie m *)+
G M (R-RM)/( (R-RM).(R-RM))^(3/2) (* Siła grawitacji newtonowskiej pochodząca od ciała o masie M *)

Sprawdzamy, czy otrzymane układy równań są równoważne.

First/@sysROT2-sysROT3//Simplify

Poniżej przykładowe rozwiązanie numeryczne w układzie korotującym.

G=1

M=10

m=1

d=1

Ω=Sqrt[G(m+M)/d^3]

rozw=NDSolve[Join[sysROT2,{X[0]==1,Y[0]==1/10
,X'[0]==0,Y'[0]==24/10}],{X,Y},{t,0,33},Method->"Automatic",WorkingPrecision->16,MaxSteps->10^6][[1]]

x=X/.rozw;

y=Y/.rozw;

ParametricPlot[{x[t],y[t]},{t,0,33},PlotRange->1.2d{{-1,1},{-1,1}},PlotStyle->Thin,PlotPoints->256,Epilog->{{Black,Point[Rm]},{PointSize[Large],Point[RM]}}]

Dalsza część wykładu będzie poszukiwaniem reguł i ograniczeń w powyżej widocznym, pozornie chaotycznym ruchu.