Mathematische Grundlagen zur umfangreichen Programmieraufgabe:

1.Teil: Diskrete Berechnung der Massenpunkte:

Aufgabe: Diskretisierung des n-Körperproblems und graphische Umsetzung. Voraussetzung: Für die diskrete Berechnung benötigen wir einen Modellraum mit folgenden Eigenschaften: -Es gelten nur newtonsche Gesetze (keine relativ.-Gesetze). -Alle Massen sind Punktförmig (minimaler Abstand und große Kräfte möglich). -Für unser Problem nur zwei-dimensional. -Energieerhaltungssatz und Wechselwirkungsgesetz sind erfüllt. Größen: n ist die Anzahl der Massenpunkte. Für den i-ten Massenpunkt bezeichne m_{i} seine Masse, x_{i} (t) seinen Ortsvektor zur Zeit t, x_{i}^{·} (t) seinen Geschwindigkeitsvektor z. Z. t, x_{i}^{··} (t) seinen Beschleunigungsvektor z. Z. t und F_{i} (t) die auf ihn wirkende gerichtete Kraft z. Z. t. Für das gesamte System bezeichne m := ∑_{i=1}^{n} m_{i} die Gesamtmasse,

F (t) := ∑_{i=1}^{n} F_{i} (t) die Kraftsumme z. Z. t und

T (t) := ∑_{i=1}^{n} \frac{m_{i}}{2} |x_{i}^{·} (t)|^{2} die gesamte kinetische Energie z. Z. t.

Außerdem bezeichne F_{ij} (t) die zwischen m_{i} und m_{j} wirkende Kraft, T die kinetische Energie Gesetze:

[N] : m_{i} × x_{i}^{··} (t) = F_{i} (t) [W W G] : F_{ji} (t) = −F_{ij} (t)

[EEG] : T_{k+1} + U_{k+1} = T_{k} + U_{k} Damit der Energieerhaltungssatz gilt, müssen folgende Voraussetzungen erfüllt sein: -Es gelten die Grundgesetze der newtonschen Mechanik. -Es existiert ein konservatives Kraftfeld. Ein konservatives Kraftfeld liegt dann vor, wenn sich die Kräfte als negative Gradienden eines Potentials Φ schreiben lassen. D. h. es ex. ein differenzierbares Φ Φ: (R^{2})^{n} ⊃ D ∋ (x_{1} (t) , ..., x_{n} (t)) ⟶ Φ (x_{1} (t) , ..., x_{n} (t)) ∈ R so daß F_{i} (t) = −grad_{i}Φ (x_{1} (t) , ..., x_{n} (t)) für alle i. Definiere grad_{i}Φ wie folgt: Φ (x_{1} (t) , ..., x_{i−1} (t) , x_{i} (t) + h, x_{i+1} (t) , ..., x_{n} (t)) = = Φ (x_{1} (t) , ..., x_{n} (t)) + grad_{i}Φ (x_{1} (t) , ..., x_{n} (t)) × h + δ (x_{1} (t) , ...x_{n} (t) , h) mit

lim_{0≠h→o} \frac{1}{|h|}δ (x_{1} (t) , ..., x_{n} (t) , h) = 0

wobei grad_{i}Φ = ( \frac{dΦ}{dx_{i,x}} ; \frac{dΦ}{dx_{i,y}} )

Konsequenz für T (t) = ∑_{i=1}^{n} \frac{m_{i}}{2} |x_{i}^{·} (t)|^{2}

T^{·} (t) = ∑_{i=1}^{n} m_{i}x_{i}^{··} (t) x_{i}^{·} (t) = ∑_{i=1}^{n} F_{i} (t) x_{i}^{·} ⇔

⇔ T^{·} (t) = − ∑_{i=1}^{n} grad_{i}Φ (x_{1} (t) , ..., x_{n} (t)) x_{i}^{·} (t) = −U^{·} (t)

für U (t) := Φ (x_{1} (t) , ..., x_{n} (t)) ,siehe def. von Φ. Weil T^{·} (t) + U^{·} (t) = 0 für alle t ⇒ T (t) + U (t) = const. für alle t ⇒ T (t) + U (t) = T (t + Δt) + U (t + Δt) E.E.G. ⇒ U (t) = Φ (x_{1} (t) , ..., x_{n} (t)) =potentielle Energie. Distanzabhängige Wechselwirkungskräfte: Sei F_{ij} (t) = f_{ij} (t) \frac{r_{ij}(t)}{|r_{ij(t)}|} mit

r_{ij} (t) := x_{j} (t) − x_{i} (t) ≠ 0 und f_{ij} (t) = |F_{ij} (t)| Annahme: Es ex. γ_{ij} (|r_{ij} (t)|) , γ_{ij} : ]0, ∞[ ⟶ R stetig diff.-bar, so daß γ_{ij} (|r_{ij} (t)|) = f_{ij} (t) für alle i,j=1,...n ; i≠j Behauptung:

Dann ex. Φ so daßF_{i} (t) := ∑_{\begin{matrix}j=1 \\ j≠i \end{matrix}}^{n} F_{ij} (t) = −grad_{i}Φ (x_{1} (t) , ..., x_{n} (t))

Beweis: Γ_{ij} sei Stammfunktion von γ_{ij} (ex. weil γ stetig diff.-bar). Mit D := {(x_{1}, ..., x_{n}) ∈ (R^{2})^{n} : x_{i} ≠ x_{j} für i ≠ j} ergibt sich Φ wie folgt:

Φ : D ∋ (x_{1}, ..., x_{n}) → \frac{1}{2} ∑_{\begin{matrix}i,j=1 \\ i≠j \end{matrix}}^{n} Γ_{ij} (|x_{j} − x_{i}|) ⇒ Beh.

Damit brauchen wir nur noch zu zeigen, daßbei gravitierenden Massen ein solches Γ existiert.

Bei gravitierenden Massen gelten folgende Kräfte: F_{ij} (t) = \frac{Gm_{i}m_{j}}{|r_{ij}(t)|^{2}} \frac{r_{ij}(t)}{|r_{ij}(t)|} = γ_{ij} (|r_{ij} (t)|) \frac{r_{ij}(t)}{|r_{ij}(t)|}

⇒ γ_{ij} (x) = \frac{Gm_{i}m_{j}}{x^{2}} für alle x>0 ⇒ Γ_{ij} (x) = −\frac{Gm_{i}m_{j}}{x} Dieses γ_{ij} erfüllt alle geforderten Voraussetzungen. Daraus ergibt sich die potentielle Energie wie folgt:

U (t) := \frac{1}{2} ∑_{\begin{matrix}i,j=1 \\ j≠i \end{matrix}}^{n} Γ_{ij} (|x_{j} (t) − x_{i} (t)|) = −\frac{G}{2} ∑_{\begin{matrix}i,j=1 \\ i≠j \end{matrix}}^{n} \frac{m_{i}m_{j}}{|r_{ij}(t)|}

Energieerhaltende Diskretisierung: Vor.: Δt möglichst kleines Intervall, t^{k} = k × Δt , Startwerte bekannt. Ansatz: Aus den Näherungen x_{i}^{k} für x_{i} (t^{k}) und v_{i}^{k} für x_{i}^{·} (t^{k}) möglichst gute Näherungen x_{i}^{k+1} und v_{i}^{k+1} zu berechnen. Hierfür werden folgende Formeln benutzt:

[i] x_{i}^{k+1} := x_{i}^{k} + dv_{i}^{k} + \frac{d^{2}}{2} \frac{1}{m_{i}} F_{i}^{k+ \frac{1}{2}}

[ii] v_{i}^{k+1} := v_{i}^{k} + d\frac{1}{m_{i}} F_{i}^{k+ \frac{1}{2}} Geht man für die Diskretisierung von einem konstanten F aus, ergeben sich die Formeln i, ii aus [N] wie folgt: x_{i}^{··} (t) = \frac{F_{i}(t)}{m_{i}} = const. ∫ x_{i}^{··} (t) dt = ∫ \frac{F_{i}(t)}{m_{i}}dt x_{i}^{·} (t) = \frac{F_{i}(t)}{m_{i}}t + c_{1} ∫ x_{i}^{·} (t) dt = ∫ \frac{F_{i}^{·}(t)}{m_{i}}t + c_{1}dt x_{i} (t) = \frac{F_{i}(t)}{m_{i}} \frac{1}{2}t^{2} + c_{1}t + c_{2} Wähle Anfangsbedingung: c_{1} = v^{k}, c_{2} = x^{k}, t = Δt = d ⇒ [i] , [ii] F_{i}^{k+ \frac{1}{2}} soll dabei als konstante Approximation von F_{i} |_{[t_{k},t_{k+1}]} so gewählt werden, daßfür konservative F_{1}, ..., F_{n} der Energiesatz gilt. Damit folgt für die kinetische Energie

T^{k} := ∑_{i=1}^{n} \frac{m_{i}}{2} |v_{i}^{k}|^{2} , k = 0, 1, 2, ....

die zeitliche Änderung aus [ii] ⇒ v_{i}^{k+1} − v_{i}^{k} = d\frac{1}{m_{i}} F_{i}^{k+ \frac{1}{2}} [∗] und

aus [i] ⇒ x_{i}^{k+1} = x_{i}^{k} + dv_{i}^{k} + \frac{d^{2}}{2} \frac{1}{m_{i}} F_{i}^{k+ \frac{1}{2}} = x_{i}^{k} + \frac{d}{2} (v_{i}^{k+1} + v_{i}^{k})

⇒ \frac{v_{i}^{k+1}+v_{i}^{k}}{2} = \frac{1}{d} (x_{i}^{k+1} − x_{i}^{k}) [∗∗]

[∗, ∗∗] ⇒ \frac{1}{2} (v_{i}^{k+1} − v_{i}^{k}) (v_{i}^{k+1} + v_{i}^{k}) = \frac{1}{2} (|v_{i}^{k+1}|^{2} − |v_{i}^{k}|^{2}) = \frac{1}{m_{i}} F_{i}^{k+ \frac{1}{2}} (x_{i}^{k+1} − x_{i}^{k})

⇒ T^{k+1} − T^{k} = ∑_{i=1}^{n} \frac{m_{i}}{2} (|v_{i}^{k+1}|^{2} − |v_{i}^{k}|^{2}) = ∑_{i=1}^{n} F_{i}^{k+ \frac{1}{2}} (x_{i}^{k+1} − x_{i}^{k})

Für die potentielle Energie folgt:

U^{k} := Φ (x_{1}^{k}, ..., x_{n}^{k}) , k = 0, 1, 2, ... falls U (t) = Φ (x_{1} (t) , ..., x_{n} (t)) Daraus ergibt sich, der E.E.S. gilt nur, falls F_{i}^{k+ \frac{1}{2}} folgende Bedingungen erfüllt:

∑_{i=1}^{n} F_{i}^{k+ \frac{1}{2}} (x_{i}^{k+1} − x_{i}^{k}) = Φ (x_{1}^{k}, ..., x_{n}^{k}) − Φ (x_{1}^{k+1}, ..., x_{n}^{k+1})

(da E.E.S. ⇔ T^{k+1} − T^{k} = U^{k} − U^{k+1} unabhängig von k) Mit dieser Bedingung ergibt sich die Wahl der F_{i}^{k+ \frac{1}{2}} für distanzabhängige Zweiteilchenwechselwirkungskräfte so: Es sei die Annahme erfüllt, daß

F_{i} (t) = ∑_{\begin{matrix}j=1 \\ j≠i \end{matrix}}^{n} F_{ij} (t) = ∑_{\begin{matrix}j=1 \\ j≠i \end{matrix}}^{n} γ_{ij} |x_{j} (t) − x_{i} (t)| t\frac{x_{j}−x_{i}}{|x_{j}−x_{i}|}

mit Stammfunktion Γ_{ij} von γ_{ij} existiert, dann gilt

Φ (x_{1} (t) , ..., x_{n} (t)) := \frac{1}{2} ∑_{\begin{matrix}i,j=1 \\ j≠i \end{matrix}}^{n} Γ_{ij} (|x_{j} (t) − x_{i} (t)|)

In unsere Bedingung für E.E.-Diskretisierung eingesetzt ergibt sich:

[E.E.G] ⇔ ∑_{i=1}^{n} F_{i}^{k+ \frac{1}{2}} (x_{i}^{k+1} − x_{i}^{k}) =

= \frac{1}{2} ∑_{\begin{matrix}i,j=1 \\ j≠i \end{matrix}}^{n} Γ_{ij} (|x_{j}^{k} − x_{i}^{k}|) − \frac{1}{2} ∑_{\begin{matrix}i,j=1 \\ j≠i \end{matrix}}^{n} Γ_{ij} (|x_{j}^{k+1} − x_{i}^{k+1}|) =

= \frac{1}{2} ∑_{\begin{matrix}i,j=1 \\ j≠i \end{matrix}}^{n} (Γ_{ij} (|r_{ij}^{k}|) − Γ_{ij} (|r_{ij}^{k+1}|)) = [1]

Setzt man

∗γ_{ij}^{k+ \frac{1}{2}} := \{ \begin{matrix}γ_{ij} (|r_{ij}^{k}|) falls |r_{ij}^{k}| = |r_{ij}^{k+1}| \\ \frac{Γ_{ij}(|r_{ij}^{k+1}|)−Γ_{ij}(|r_{ij}^{k}|)}{|r_{ij}^{k+1}|−|r_{ij}^{k}|} & sonst \end{matrix} \}

dann wird

[1] = \frac{1}{2} ∑_{\begin{matrix}i,j=1 \\ j≠i \end{matrix}}^{n} (−^{∗}γ_{ij}^{k+ \frac{1}{2}} \frac{(|r_{ij}^{k+1}|−|r_{ij}^{k}|)(|r_{ij}^{k+1}|+|r_{ij}^{k}|)}{|r_{ij}^{k+1}|+|r_{ij}^{k}|})

unter Verwendung der binomischen Eigenschaft (|a| − |b|) (|a| + |b|) = |a|^{2} − |b|^{2} = (a − b) (a + b) kann man die Beträge weglassen.

[1] = \frac{1}{2} ∑_{\begin{matrix}i,j=1 \\ j≠i \end{matrix}}^{n} (−^{∗}γ_{ij}^{k+ \frac{1}{2}} \frac{(r_{ij}^{k+1}−r_{ij}^{k})(r_{ij}^{k+1}+r_{ij}^{k})}{|r_{ij}^{k+1}|+|r_{ij}^{k}|})

[r_{ij}^{k+1} − r_{ij}^{k} = (x_{j}^{k+1} − x_{j}^{k}) − (x_{i}^{k+1} − x_{i}^{k}) ⇒]

[1] = \frac{1}{2} ∑_{\begin{matrix}i,j=1 \\ j≠i \end{matrix}}^{n} (−^{∗}γ_{ij}^{k+ \frac{1}{2}} (x_{j}^{k+1} − x_{j}^{k}) \frac{r_{ij}^{k+1}+r_{ij}^{k}}{|r_{ij}^{k+1}|+|r_{ij}^{k}|} +^{∗} γ_{ij}^{k+ \frac{1}{2}} (x_{i}^{k+1} − x_{i}^{k}) \frac{r_{ij}^{k+1}+r_{ij}^{k}}{|r_{ij}^{k+1}|+|r_{ij}^{k}|}) ⇒

[1] = \frac{1}{2} ∑_{\begin{matrix}i,j=1 \\ j≠i \end{matrix}}^{n} (−^{∗}γ_{ij}^{k+ \frac{1}{2}} (x_{j}^{k+1} − x_{j}^{k}) \frac{r_{ij}^{k+1}+r_{ij}^{k}}{|r_{ij}^{k+1}|+|r_{ij}^{k}|})+\frac{1}{2} ∑_{\begin{matrix}i,j=1 \\ j≠i \end{matrix}}^{n} (+^{∗}γ_{ij}^{k+ \frac{1}{2}} (x_{i}^{k+1} − x_{i}^{k}) \frac{r_{ij}^{k+1}+r_{ij}^{k}}{|r_{ij}^{k+1}|+|r_{ij}^{k}|})

⇒ [1] = ∑_{\begin{matrix}i,j=1 \\ j≠i \end{matrix}}^{n} (+^{∗}γ_{ij}^{k+ \frac{1}{2}} (x_{i}^{k+1} − x_{i}^{k}) \frac{r_{ij}^{k+1}+r_{ij}^{k}}{|r_{ij}^{k+1}|+|r_{ij}^{k}|})

[folgt aus einfacher Umordnung der Reihen mit r_{ij}^{k} = −r_{ji}^{k} und ^{∗}γ_{ij}^{k+ \frac{1}{2}} =^{∗} γ_{ji}^{k+ \frac{1}{2}} ] ⇒der E.E.S. ist damit äquivalent zu

∑_{i=1}^{n} F_{i}^{k+ \frac{1}{2}} (x_{i}^{k+1} − x_{i}^{k}) = ∑_{i=1}^{n} ( ∑_{\begin{matrix}j=1 \\ j≠i \end{matrix}}^{n} +^{∗}γ_{ij}^{k+ \frac{1}{2}} \frac{r_{ij}^{k+1}+r_{ij}^{k}}{|r_{ij}^{k+1}|+|r_{ij}^{k}|}) (x_{i}^{k} − x_{i}^{k+1})

Setzt man nun F_{ij}^{k+ \frac{1}{2}} :=^{∗} γ_{ij}^{k+ \frac{1}{2}} \frac{r_{ij}^{k+1}+r_{ij}^{k}}{|r_{ij}^{k+1}|+|r_{ij}^{k}|} und

F_{i}^{k+ \frac{1}{2}} := ∑_{\begin{matrix}j=1 \\ j≠i \end{matrix}}^{n} F_{ij}^{k+ \frac{1}{2}}

dann ist unsere Energieerhaltung erfüllt. Daraus ergibt sich für unsere Berechnung aus Γ_{ij} (x) = −\frac{Gm_{i}m_{j}}{x} für x ≠ y :

∗γ_{ij}^{k+ \frac{1}{2}} = \frac{Γ_{ij}(x)−Γ_{ij}(y)}{x−y} = −Gm_{i}m_{j} ( \frac{\frac{1}{x} − \frac{1}{y}}{x−y} ) = Gm_{i}m_{j} \frac{1}{xy} = Gm_{i}m_{j} _{\frac{1}{|r_{ij}^{k+1}||r_{ij}^{k}|}}

und F_{ij}^{k+ \frac{1}{2}} = Gm_{i}m_{j} _{\frac{1}{|r_{ij}^{k+1}||r_{ij}^{k}|}} \frac{r_{ij}^{k+1}+r_{ij}^{k}}{|r_{ij}^{k+1}|+|r_{ij}^{k}|}

Aus diesen F_{ij}^{k+ \frac{1}{2}} ergeben sich dann letztlich folgende Formeln zur Berechnung der neuen Geschwindigkeiten und Positionen.

x_{i}^{k+1} := x_{i}^{k} + dv_{i}^{k} + \frac{d^{2}}{2} ∑_{\begin{matrix}j=1 \\ j≠i \end{matrix}}^{n} \frac{m_{j}G}{|r_{ij}^{k+1}| |r_{ij}^{k}|} \frac{r_{ij}^{k+1} + r_{ij}^{k}}{|r_{ij}^{k+1}| + |r_{ij}^{k}|}

v_{i}^{k+1} := v_{i}^{k} + d ∑_{\begin{matrix}j=1 \\ j≠i \end{matrix}}^{n} \frac{m_{j}G}{|r_{ij}^{k+1}| |r_{ij}^{k}|} \frac{r_{ij}^{k+1} + r_{ij}^{k}}{|r_{ij}^{k+1}| + |r_{ij}^{k}|}

2.Teil: Schrittweitensteuerung:

Aufgabe: Da die Schrittweite sehr stark die Genauigkeit der Berechnung beeinflußt, ist es sinnvoll, die jeweilige Schrittweite von der gegebenen Situatuon bzw. Genauigkeit abhängig zu machen. Bei konstanten Werten entsteht das Problem einer zu kleinen Schrittweite, wenn die Bewegungen der Planeten fast gradlinig gleichförmig sind und einer zu großen Schrittweite (d.h. auch zu große Ungenauigkeit) bei nicht gleichförmigen Bewegungen der Massen. Aus diesem Grund wollen wir eine Schrittweitensteuerung einbauen und bedienen uns dabei des Runge-Kutta-Verfahrens [vgl. J.Stoer / Numerische Mathematik]. Für Interessierte an den mathematischen Grundlagen der Steuerung verweisen wir auf das Buch und den Auszug im Anhang. Umsetzung: Die Genauigkeit ϵ wird im Programm eingegeben und ist bekannt. 1: Schrittweite Δt aus vorherigem Durchlauf bekannt. 2: Aus x (t) wird x (t + 1) auf zwei Wegen berechnet. Einmal xa (t + 1) mit Δt und einmal xb (t + 1) in zwei Durchgängen mit jeweils \frac{Δt}{2} als Schrittweite. 3: Die Werte werden verglichen. D.h. es wird die Summe der absoluten Differenzen der Sternpositionen als b errechnet. 4: Berechne (i) hd := \sqrt{\frac{2b}{nϵ}} (n Anzahl der Massen) (ii) h := 2\frac{Δt}{hd} und (iii) Δt = 2h . 5: Wenn hd ≤ 2 dann setze x (t + 1) = xb und berechne x (t + 2) mit dem neuen, größeren Δt. Sonst rechne erneut mit altem x (t) und dem verkleinertem Δt ab (1:). Dies wird wiederholt, bis die Bedingung erfüllt wird.