Gelöste Aufgaben/T3BP

Aufgabenstellung

Wir untersuchen das „Three-Body-Problem“(vgl. Wikipedia) numerisch. Dabei sollen die Bahnen von drei Körper mit den Punktmassen m₁, m₂, m₃ in Wechselwirkung miteinander berechnet werden.

"Die Drei Sonnen"

Gesucht ist die Lösung des Anfangswertproblems für verschiedene Anfangswerte (Orte und Geschwindigkeiten) und Massen m_i der Körper.

Lösung mit Matlab®

Die Lösung mit Matlab erfordert keine großen algebraischen Vorbereitungen - für die wir sonst Maxima einsetzten würden. Als äußere Kräfte treten nur die Feldkräfte der Gravitation zwischen den Körpern auf - die wir einfach hinschreiben und in Matlab implementieren können.

Declarations

Dafür verwenden wir das script T3BP.m, das die Funktionen

• preprocess.m,
• ode45.m (solve) und
• postprocess.m

aufruft. Die Klasse

• planets.m

verwenden wir nur, um die Systemparameter "sys" bequem ansprechen zu können. Zur Lösung des Anfangswertproblems rufen wir eine der ODE-Funktionen aus Matlab auf, die wir mit den Gleichgewichtsbedingungen unserer Funktion t3bpdydt.m füttern.

T3BP
   |- preprocess
          |- class sys=planets(data)
   |- ode45
          |- t3bpdydt
   |- postprocess

Equilibrium Conditions

Die Gravitationskräfte zwischen den zwei Körpern i und j erfasst

$F_{i,j}=G\cdot {\frac {\displaystyle m_{i}\cdot m_{j}}{\displaystyle r_{i,j}^{2}}}$

mit der Gravitationskonstanten

$G=6.674\cdot 10^{-11}{\frac {\displaystyle m^{3}}{\displaystyle kg\cdot s^{2}}}.$

Dabei gehen die Wirkungslinien der Kräfte durch die Massenmittelpunkte der Körper, die Körper ziehen sich gegenseitig an. Damit ist

${\vec {F}}_{i,j}=F_{i,j}\cdot {\vec {e}}_{i,j}$,

wobei ${\vec {e}}_{i,j}$ der Einheitsvektor der Anziehungskraft - in diesem Fall von i nach j - ist.

Die Bewegungsgleichungen schreiben wir mit $F$ in Vektorschreibweise für Körper i als

$m_{i}{\dot {\ddot {\vec {u}}}}_{i}=\sum _{\ell =j,k}G\cdot \underbrace {\frac {\displaystyle m_{i}\cdot m_{\ell }}{\displaystyle r_{i,\ell }^{2}}} _{=F_{i,j}}\cdot {\vec {e}}_{i,\ell }{\text{ mit }}{\vec {u}}_{i}=\left({\vec {e}}_{x},{\vec {e}}_{y},{\vec {e}}_{z}\right)\cdot \underbrace {\left({\begin{array}{l}u_{x}\\u_{y}\\u_{z}\end{array}}\right)} _{\displaystyle ={\underline {u}}_{i}}$.

Dabei ist z.B. i=1 und $\ell$=2,3.

Die drei Bewegungsgleichungen in den drei räumlichen Koordinaten u_i formulieren wir in dimensionslosen Koordinaten. Dafür brauchen wir drei unabhängige Referenzgrößen, hier wählen wir

${\begin{array}{lcll}M&=&m_{1}&{\text{ Referenz-Masse}}\\L&=&1.610^{11}{\text{km}}&{\text{ Referenz-Länge: der Durchmesser unseres Sonnensystems und}}\\F&=&G\cdot {\frac {\displaystyle m_{1}\cdot m_{1}}{\displaystyle L^{2}}}&{\text{ Referenz-Kraft}}\end{array}}.$

Uns fehlt noch die Referent-Zeit T, die wir aus

$F={\frac {\displaystyle m_{1}\cdot L}{\displaystyle T^{2}}}{\text{ zu }}T={\sqrt {\frac {\displaystyle m_{1}\cdot L}{\displaystyle F}}}$

erhalten. Damit ist $T\approx 175Jahre$.

Damit können wir schreiben:

${\begin{array}{lcll}t&=&\tau \cdot T&{\text{ mit der dimensionslosen Zeit }}\tau \\m_{i}&=&\theta _{i}\cdot m_{1}&{\text{ mit der dimensionslosen Masse }}\theta _{i}\\{\underline {u}}_{i}(t)&=&L\cdot {\underline {U}}_{i}(\tau )&{\text{ mit den dimensionslose Koordinaten }}U_{i}{\text{ und}}\\r_{i,j}&=&L\cdot \varrho _{i,j}&{\text{ mit dem dimensionslosen Abstand zweier Körper }}\varrho _{i,j}\end{array}}.$

Einsetzen und Kürzen liefert uns dann die dimensionslosen Bewegungsgleichungen

${\ddot {\underline {U}}}_{i}=\sum _{\ell =j,k}{\frac {\displaystyle \theta _{\ell }}{\displaystyle \varrho _{i,\ell }^{2}}}\cdot {\underline {e}}_{i,\ell }$.

Diese Gleichgewichtsbeziehungen implementieren wir in der Funktion t3bpdydt.m, die vom Löser ODE45 aufgerufen wird.

function dydt = t3bpdydt(t,y,sys)
% implementation of ode
% sys holds system parameters
% get coordinates
for body = 1:3
    u(body,1:3) = transpose(y(3*(body-1)+1:3*body,1));
end

r = zeros(3,3,3);
e = zeros(3,3,3);

% compose vectors, unit vectors and distances
r(1,2,:) = u(2,:)-u(1,:); % this from m[1] to m[2]
   :
   :
   :
% mass 3
dydt(9+7:9+9,1)=sys.theta(1)/R(3,1)*e(3,1,:)+sys.theta(2)/R(3,2)*e(3,2,:);

%% velocities
dydt(1:9,1) = y(10:18,1);

%%
waitbar(t / sys.tEnd);

end

Solving

Die Lösung des Anfangswertproblems übertragen wir der Matlab-Routine ODE45. Dabei organisieren wir die Zustandsgrößen so, dass in y0 zunächst die 9 Anfangsverschiebungen und dann die 9 Anfangsgeschwindigkeiten stehen.

Wir lösen dann

${\dot {\underline {q}}}(\tau )={\underline {f}}({\underline {q}}(\tau ))$

mit

$\underbrace {{\underline {q}}(0)} _{={\text{y0}}}=\left({\begin{array}{c}U_{1,x}(0)\\U_{1,y}(0)\\U_{1,z}(0)\\U_{2,x}(0)\\\vdots \\V_{3,y}(0)\\V_{3,z}(0)\end{array}}\right).$

Die $V_{i,\xi }$ sind Komponentenweise die Geschwindigkeiten, also z.B.

$V_{3,y}(\tau )={\frac {\displaystyle d}{\displaystyle d\tau }}U_{3,y}(\tau )$.

%% solve IVP
% set-up numerical paramters
tspan = 0:sys.tEnd/1000:sys.tEnd;
u0 = [transpose(sys.u0(1,:));transpose(sys.u0(2,:));transpose(sys.u0(3,:))]/sys.L;
v0 = [transpose(sys.v0(1,:));transpose(sys.v0(2,:));transpose(sys.v0(3,:))]*sys.T/sys.L;
% initial values from sys structured for ODE45
y0    = [u0;v0];
h = waitbar(0,'solving ODE - please wait...');
%
[t,y] = ode45(@(t,y)t3bpdydt(t,y,sys),tspan,y0);
close(h);

Matlab^©-files

Folder Structure

Die Datei-Struktur des Projekts zeigt das Bild links mit dem Skript T3BP.m. Classes und Functions sind in den jeweiligen Ordnern. Die Excel-Datei hält alle System-Parameter.

Den komplette Quellcode zu diesem Programm können Sie über dieses ZIP-File rechts herunterladen.

archive

download compressed archive →

Post-Processing

Trajektoren der Körper

Uns interessieren natürlich die Trajektorien der Körper im Raum - dabei zeichnen wir die Körper-Durchmesser entsprechend Ihrer Masse. Die Darstellung ist nicht maßstäblich.

Zur Kontrolle können wir prüfen, ob die ODE-Routine im Rahmen der numerischen Genauigkeit die Erhaltungssätze der Mechanik erfüllt. Das machen wir hier am Beispiel des Impulssatzes: die Summe der Bewegungsgrößen jeweils in die drei Raumrichtungen muss gleich bleiben:

Bewegungsgrößen $\sum _{i=1}^{3}\theta _{i}\cdot I_{i,x0}$

.

Und das tut sie offensichtlich auch - im Rahmen der numerischen Genauigkeit.

Und schließlich schauen wir uns die Bewegung der drei Körper in einer Animation über der Zeit an - hier über einen Zeitraum 4375 Jahren.

Video laden

Lokale Datei

Lokale Datei sammelt möglicherweise persönliche Daten.

FortfahrenAusblenden

Animation der Bewegung

nothing here

---

Literature

• The Three-Body-Problem