Gelöste Aufgaben/T3BP

Aufgabenstellung

Sie 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 Maxima einsetzten würden. Als äußere Kräfte treten nur die Feldkräfte der Gravitation zwischen den Körpern auf - die wir einfach einschreiben und in Matlab implementieren können.

Declarations

Wir implementieren die Lösung in Matlab. 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.

📑 toggle code

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

Equilibrium Conditions

Die Gravitationskräfte zwischen den Körpern i und j erfassen wir mit

$F=G\cdot \frac{m_i\cdot m_j}{r_{i,j}^2}$

mit der Gravitationskonstanten

$G=6.674\cdot 10^{-11}\frac{m^3}{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

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

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

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

$m_i{\dot{\ddot{\overrightarrow{u}}}}_i=\sum _{\ell =j,k}G\cdot \frac{m_i\cdot m_{\ell }}{r_{i,\ell }^2}{\dot{\overrightarrow{e}}}_{i,\ell }\text{ mit }{\overrightarrow{u}}_i=({\overrightarrow{e}}_x,{\overrightarrow{e}}_y,{\overrightarrow{e}}_z)\cdot \underset{={\underset{\_}{u}}_i}{\underbrace{\left(\begin{array}{l}{u}_x\\ u_y\\ u_z\end{array}\right)}}$.

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{m_1\cdot m_1}{L^2}& \text{ Referenz-Kraft}\end{array}.$

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

$F=\frac{m_1\cdot L}{T^2}\text{ zu }T=\sqrt{\frac{m_1\cdot L}{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\\ {\underset{\_}{u}}_i(t)& =& L\cdot {\underset{\_}{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{\underset{\_}{U}}}_i=\sum _{\ell =j,k}\frac{{\theta }_{\ell }}{{\varrho }_{i,\ell }^2}{\dot{\underset{\_}{e}}}_{i,\ell }$.

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

📑 toggle code

function dydt = t3bpdydt(t,y,sys)
% implementation of ode
% params hold 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{\underset{\_}{q}}(\tau )=\underset{\_}{f}(\underset{\_}{q}(\tau ))$

mit

$\underset{\_}{q}(0)=\left(\begin{array}{c}{U}_{1,x}(0)\\ U_{1,y}(0)\\ U_{1,z}(0)\\ U_{2,x}(0)\\ ⋮\\ V_{3,y}(0)\\ V_{3,z}(0)\end{array}\right).$

z.B. mit

Fehler beim Parsen (Syntaxfehler): {\displaystyle V_{3,y}(\tau) = \frac{\displaystyle d}{\displaystyle d\tau} U_{3,y}(\tau). |code= <syntaxhighlight lang="matlab" line start=1> %% 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); </syntaxhighlight> }} {{MyZipfileBlock |title=Matlab^©-files |text =[[Datei:matlab-folder-structure.png|links|mini|Folder Structure]]Die Datei-Struktur zeit das Skript <nowiki>T3BP.m</nowiki>. 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. |file =T3BP.zip}} <!-------------------------------------------------------------------------------->{{MyCodeBlock|title=Post-Processing |text=[[Datei:T3BP-21.png|550px|right|mini|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 Impulssystes: die Summe der Bewegungsgrößen jeweils in die drei Raumrichtungen muss gleich bleiben: [[Datei:T3BP-22.png|150px|left|mini|Bewegungsgrößen <math>\sum_{i=1}^3 \theta_i \cdot I_{i,x0}} ]].

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 abspielen

Lokale Datei

Lokale Datei sammelt möglicherweise persönliche Daten.

FortfahrenAusblenden

Animation der Bewegung

📑 toggle code

1+1

---

Links

• The Three-Body-Problem

Literature

• ...