kleiner planetarer Simulator

Bernhard

Registriertes Mitglied
Bei einer Schrittweite von 10000 Sekunden ergibt sich:
y_n + 0.5 * h * k1 = (-7.48e10, 1.0525e8, 0.0, 2.9654e7, 2.105e4, 0.0, 7.48e10, -1.0525e8, 0.0, -2.9654e7, -2.105e4, 0.0)
Muss demnach zu


Bei einer Schrittweite von 10000 Sekunden ergibt sich:
y_n + 0.5 * h * k1 = (-7.48e10, 1.0525e8, 0.0, 2.9654e1, 2.105e4, 0.0, 7.48e10, -1.0525e8, 0.0, -2.9654e1, -2.105e4, 0.0)


korrigiert werden.
 

Kibo

Registriertes Mitglied
Wenn ich deinen Vektor richtig verstanden habe, dann ist bei dir

Die Position des Prognosepunkts P1 von Sonne1 (-7.48e10/ 1.0525e8/ 0.0)
und die "Steigung" m1 dann (2.9654e1, 2.105e4, 0.0)

Die von meinem Programm errechneten Werte sind ja
für P1 (-74799998516,438034.../10525000/0)
und für m1 (2,967124.../21050/0)

Es unterscheiden sich also der x-Wert von P1 und m1. Sehr wundert mich ersterer, da du selber schreibst

Die Beschleunigung als Vektor zeigt am Anfang also in Richtung der positiven x-Achse.

Also warum prognostiziert dein P1 diese Änderung nicht?

Wie wir übereinstimmend feststellten ist die durch die Gravitation verursachte Beschleunigung (rund) 5.93e-3 m/s^2, diese wirkt nur in Richtung x Achse. Da wir das für die halbe gewählte Schrittweite (500 s) betrachten kommen wir beide auf diese nahezu gleichen Werte, warum sie allerdings nicht 100%tig gleich sind, wüsste Ich schon gerne. Eine Kontrollrechnung in einer Exceltabelle, bestätigt die Werte meines Programms.

Nichtsdestotrotz muss mein Runge Kutta falsch sein, da bei einem Probelauf sich die Orbits der Sonnen schon nach kurzer Zeit immer weiter absenken sodass es zu Kollision kommt (Energieerhaltung???).
Euler hingegen, liefert selbst bei 800000 Sekunden Schrittweite noch nahezu perfekt kreisförmige Orbits, nach simulierten 1000 Jahren ist keine Destabilisierung erkennbar.

mfg
 
Zuletzt bearbeitet:

Bernhard

Registriertes Mitglied
Also warum prognostiziert dein P1 diese Änderung nicht?
Die Gravitationskraft ist eine Beschleunigung und diese ändert ihrerseits im ersten Schritt nur die Geschwindkeit. Der Prognosepunkt enthält eben auch die Geschwindigkeit des Massepunktes 1 und die x-Komponente dieser Geschwindigkeit wird dann ja auch auf den Wert 29.7 m/s gesetzt.

warum sie allerdings nicht 100%tig gleich sind, wüsste Ich schon gerne.
Ich habe 10000 (zehntausend) Sekunden als Schrittweite verwendet. Wenn ich 1000 Sekunden verwende bekomme ich genau die gleichen Werte. Der Fehler liegt vermutlich wo anders. Deswegen rechne ich einfach mal weiter:

k2 = (2.9654e1, 2.105e4, 0.0, 5.9308e-3, -8.3452e-6, 0.0, -2.9654e1, -2.105e4, 0.0, -5.9308e-3, 8.3452e-6, 0.0)

sowie

y_n + 0.5 * h * k2 = (-7.48e10, 1.0525e8, 0.0, 2.9654e1, 2.105e4, 0.0, 7.48e10, -1.0525e8, 0.0, -2.9654e1, -2.105e4, 0.0)

und damit k3 gleich k2 in den hier aufgeschriebenen Dezimalstellen.
 
Zuletzt bearbeitet:

Bernhard

Registriertes Mitglied
Hallo Kibo,

Da die dortige Rechnung eine Änderung der Koordinaten impliziert, ist das dann dort falsch dargestellt?

das ist, so weit ich es verstehe, ein Mißverständnis. In der englischen Wikipedia wird anstelle von m1 bis m4 nur die Bezeichnung k1 bis k4 verwendet. Der Algorithmus ist ansonsten der selbe. Die Koordinaten sind immer euklidisch. Man muss nur beachten, dass man bei einem Zweikörperproblem die zwei Orte x_1 und x_2 und zusätzlich noch die zwei Geschwindigkeiten v_1 und v_2 praktisch als unabhängige Variablen betrachten muss. x_1, x_2, v_1 und v_2 sind dabei die bekannten, dreidimensionalen Vektoren im euklidischen Raum, die den Bewegungszustand der beiden Körper charakterisieren.

m1 bis m4 hat also jeweils zwölf Komponenten, bzw. Einträge. Bei den Prognosepunkten werden ebenfalls für alle vier Vektoren neue Werte berechnet und zwar wieder als zwölf reelle Zahlen (4 * 3 = 12).
MfG
 

Kibo

Registriertes Mitglied
Nach einigem Kopfzerbrechen, gibt's ein Paar Fortschritte. Die Bahnen sind jetzt auch im RK Kreisrund wie sie sein sollten. Es existiert abe rnoch mindestens 1 weiterer Programmierfehler, da das System immernoch mit der Zeit seine Bahnenergie verliert (Todesspirale).
 

Bernhard

Registriertes Mitglied
Ich vermute, dass immer noch ein Fehler in der Verwaltung, bzw. softwareseitigen Trennung der verschiedenen Variablen vorliegt. Bedenke, dass die Gravitationskraft nur die aktuellen Geschwindigkeiten ändert. Die Positionen (Orte) werden dagegen nur über die Geschwindigkeiten beeinflusst.
 

Kibo

Registriertes Mitglied
Die 4 Sonnen, kreisen auch wieder wie gewohnt, es gab Inkompatibilitäten mit meiner Originalschleife und newton nr3. 1.7 rückt in greifbare nähe, jetzt muss ich noch den Spiralenfehler beseitigen.

mfg
 

Bernhard

Registriertes Mitglied
Hallo Kibo,

anbei die letzten Zahlen für den ersten Schritt:

y + h * k_3 = (-7.48e10, 2.105e8, 0.0, 5.9308e1, 2.105e4, 0.0, 7.48e10, -2.105e8, 0.0, -5.9308e1, -2.105e4, 0.0)

k_4 = (5.9308e1, 2.105e4, 0.0, 5.9308e-3, -1.669e-5, 0.0, -5.9308e1, -2.105e4, 0.0, -5.9308e-3, 1.669e-5, 0.0)

Der neue Gesamtvektor (Ort + Geschwindigkeit für beide Sonnen) ergibt sich dann zu

y_1 = (-7.48e10, 2.105e4, 0.0, 5.9308e-3, 2.105e4, 0.0, 7.48e10, -2.105e4, 0.0, -5.9308e-3, -2.105e4, 0.0)

bei einer Schrittweite von 10.000 Sekunden bekommen die beiden y-Koordinaten eine Korrektur und die beiden Geschwindigkeitsvektoren bekommen bei der x-Komponente eine Korrektur. Die vier Korrekturen bestehen dabei aus zwei Werten, die mit beiden Vorzeichen eingehen.
 

Kibo

Registriertes Mitglied
Also Bernhard, wenn ich das jetzt richtig verstanden habe, wird in einem Schritt bei dir folgendes gemacht:

Die Koordinaten der Körper werden aufgrund der Geschwindigkeit angepasst

DANACH

werden die Beschleunigungen ausgerechnet und die Geschwindigkeiten angepasst.

Sehe ich das so richtig? Wenn ja, warum nicht anders herum?
 

Bernhard

Registriertes Mitglied
Sehe ich das so richtig?
So in etwa. Ich erkläre es aber gerne auch etwas genauer:

Die Gleichung F = m * a reicht auch zusammen mit dem newtonschen Gravitationsgesetz nicht aus, um Runge-Kutta zu implementieren, denn Runge-Kutta kann nur erste Ableitungen nähern. Da die Beschleunigung nun aber ungünstigerweise die zweite Ableitung des Ortes ist kann man nicht ausschließlich mit der Beschleunigung arbeiten.

Man verwendet vielmehr den Umweg über die Geschwindigkeit und das wie folgt: Die erste Ableitung des Ortes ist die Geschwindigkeit. Ferner gilt noch: Die erste Ableitung der Geschwindigkeit ist die Beschleunigung und die kann/soll/muss man aus dem newtonschen Gravitationsgesetz ausrechnen.

Dann denkt man sich den Ort und die Geschwindigkeit als eine große Variable (mit mehreren Komponenten) und die erste Ableitung dieser Variable ist jetzt bekannt, wie im zweiten Absatz eben gezeigt. Man muss dann nur noch zwischen der Ableitung und dem Wert der großen Variablen peinlich genau unterscheiden und ist vom Verständnis her fertig.

Zur Wiederholung:
a) In der zugrundeliegenden Differentialgleichung taucht die Geschwindigkeit einmal als Wert auf und gibt dort die erste Ableitung des Ortes an und
b) Die erste Ableitung der Geschwindigkeit wird gemäß F/m berechnet, wobei F die newtonsche Gravitationskraft auf den betreffenden Körper ist.
MfG
 
Zuletzt bearbeitet:

Kibo

Registriertes Mitglied
Hmm,

also ich fass meine jetzige Lösung vielleicht am besten nochmal in Pseudocode zusammen:

Schleife bis alle Planeten durch sind
{
sichern des momentanen Geschwindigkeitsvektors in PlanetA.h & PlanetA.h

m1a[x,y,z] = gravitation(PlanetA koordinaten, PlanetB)
m1v = h + m1a*schrittweite Kommentar: Das entspricht nicht deiner Beschreibung, aber wie anders?
P1 = PlanetA koordinaten + m1v* Schrittweite*0,5
m2a[x,y,z] = gravitation(P1, PlanetB)
m2v = h + m2a*schrittweite
P2 = PlanetA koordinaten + m2v* Schrittweite*0,5
m3a[x,y,z] = gravitation(P2, PlanetB)
m3v = h + m3a*schrittweite
P3 = PlanetA koordinaten + m3v* Schrittweite
m3a[x,y,z] = gravitation(P3, PlanetB)
m3v = h + m3a*schrittweite

PlanetA.v:=( m1v + 2 * m2v + 2 * m3v + m4v)/6

PlanetA.a:= PlanetA.v - h Kommentar: Rückrechnen für Werte für PlanetB (Newton Nr3)
PlanetB.v:= -1 * PlanetA.a * PlanetA.m /PlanetB.m

}

gravitation(PlanetA koordinaten, PlanetB)
{
differenzvektor = koordninaten PlanetA - PlanetB
abstand = wurzel(Differenzvektor ²)
einheitsvektor für richtung n = differenzverktor/ abstand

Beschleunigung a = PlanetB.gm/abstand² Kommentar: PlanetB.gm = G*Masse bei Programmstart berechnet
grav Beschleunigungsvektor g = a * n
rückgabe g, abstand
}

Schleife für jeden Planeten
{
Planetkoordinaten:=Planetenkoordinaten + Planet.v*Schrittweite
}

Führt wie gesagt zu einer "Todesspirale"

MfG
 

Bernhard

Registriertes Mitglied
sichern des momentanen Geschwindigkeitsvektors in PlanetA.h & PlanetA.h
Das ist ungünstig weil verwirrend. Warum nicht PlanetA.v für die Geschwindigkeit von PlanetA und PlanetA.x für den aktuellen Ort von PlanetA? Für die Geschwindigkeit nimmt man normalerweise immer das kleine v wegen englisch velocity.

m1a[x,y,z] = gravitation(PlanetA koordinaten, PlanetB)
muss heißen
m1a[x,y,z] = gravitation(PlanetB koordinaten, PlanetA)

Planet B zieht Planeten A in Richtung hin zu B. Die Kraft auf Planet A muss also von A nach B zeigen. Da haben wir also schon mal einen Vorzeichenfehler.

m1v = h + m1a*schrittweite Kommentar: Das entspricht nicht deiner Beschreibung, aber wie anders?
muss heißen:

m1v = aktuelle Geschwindigkeit des Planeten, also beispielsweise

m1v = PlanetA.v

Die Schrittweite h wird ferner nirgends irgendwo dazu addiert, weil die Zeit nirgends explizit vorkommt. h tritt dagegen überall nur multiplikativ auf. Jedes "h + ..." ist also falsch und muss korrigiert werden.

EDIT: Bei jedem Schleifendurchlauf berechnet man aus den aktuellen Positionen und Geschwindigkeiten der Planeten und der Schrittweite die neuen Positionen und Geschwindigkeiten der Planeten. Die Zeit kann man völlig seperat berechnen oder vorher festlegen. Man kann beispielsweise ein t0, ein t_stop und eine Anzahl von Schritten N vorgeben. Dann berechnet man h = (t_stop - t0) / N. Die Werte von t0 und t_stop gehen in die Rechnung dann gar nicht mit ein. Man muss nur die Schleife genau N mal durchlaufen und darin das eben berechnete h als Schrittweite verwenden.
MfG
 
Zuletzt bearbeitet:

Kibo

Registriertes Mitglied
Guten Morgen Bernhard,

h ist die Hilfsvariable und speichert nur das alte v für die Trägheit. Die Schrittweite hab ich hier konkret mit schrittweite bezeichnet. Der letzte post von mir ist Pseudocode, und sollte eigentlich nur meine prinzipielle Vorgehensweise verständlich darstellen, ohne auf die Eigenheiten der verwendeten Programmiersprache einzugehen.
Da meine Planeten alle von einander angezogen und nicht abgestoßen werden, geh ich mal davon aus das der von dir gefundene Vorzeichenfehler auf meiner Unfähigkeit beruht meinen Code verständlich auszudrücken.:eek:

OK, du hast erst mal nur meine unverständlichen Variablennamen bemängelt, du siehst da also keinen logischen Fehler im Code an sich?

Ich lade mal das aktuell kompilierte Programm hoch (1.7alpha zu den Verbesserungen später mehr), damit du dir mal den Fehler ansehen kannst, bei 80000 Schrittweite sieht man es ganz gut. Sieht so ein bisschen nach Rundungsfehler aus, oder kann das hier ran liegen?:
Die Gleichung F = m * a reicht auch zusammen mit dem newtonschen Gravitationsgesetz nicht aus, um Runge-Kutta zu implementieren, denn Runge-Kutta kann nur erste Ableitungen nähern. Da die Beschleunigung nun aber ungünstigerweise die zweite Ableitung des Ortes ist kann man nicht ausschließlich mit der Beschleunigung arbeiten.

hier nochmal die Links:

Rapidshare

Skydrive

mfg
 
Zuletzt bearbeitet:

Bernhard

Registriertes Mitglied
Guten Morgen Kibo,

sinnvolle Bezeichnungen für Variablennamen sind extrem wichtig, damit man den Code auch nach längeren Pausen noch lesen kann. Ich schlage deswegen vor, wir suchen erst mal gemeinsam ein paar gute Bezeichnugen, um sinnvoll diskutieren zu können.

Das m1a und m1v gefällt mir schon mal gut. Es sollte für jeden Planeten exisitieren. Genauso m2a, m2v, m3a, m3v, m4a, m4v. Ferner benötigst Du für jeden Planeten einen Prognosepunkt und eine Prognosegeschwindigkeit.

Zusätzlich brauchst Du eine Funktion, welche die Differentialgleichung repräsentiert. SRMeister hatte diese Funktion in seinem Code "Evaluate" genannt. Diese Funktion erhält als Parameter die Schrittweite, den Prognosepunkt und die Prognosegeschwindigkeit und berechnet daraus die Werte für mia und miv, mit i = 1,2,3,4. Bei dieser Funktion helfe ich Dir dann gerne, wenn wir uns bei den benötigten Variablen geeinigt haben.
MfG
 

Kibo

Registriertes Mitglied
Ok, ich fasse nochmal meine jetzigen Variablen zusammen:

%Planeten%.h(x,y oder z) hilfsvariable, enthält
%Planeten%.x, y oder z Koordinaten
%Planeten%.v(x,y oder z) Geschwindigkeit
%Planeten%.g(x,y oder z) Beschleunigung

Prognosevariablen
%Planeten%.m(1, 2, 3 oder 4)v(x,y oder z)
%Planeten%.m(1, 2, 3 oder 4)a(x,y oder z)
%Planeten%.P(1, 2 oder 3)(x,y oder z)

in der Funktion (auch alternativ als Klasse probiert) zur Gravitationsberechnung
%Planeten%.(x, y oder z)d
%Planeten%..(x, y oder z)n
%Planeten%.abstand
%Planeten%.a (ohne richtung)

Das müsste alles sein, bin gerade auf Arbeit, da ist das etwas schwierig:)

%Planeten% ist eine dynamische Variable und ändert sich innerhalb der schleife entsprechend in den Namen des Aktuellen Planeten. %PlanetenB% macht das selbe für den Gegnerplaneten. Wie erwähnt nach Newton Nr 3 die Beschleunigung für %PlanetenB% gleich mit gesetzt. Es werden in einem Durchlauf auch alle Planeten genau einmal berechnet, das hab ich kontrolliert um Fehler in der Schleife auszuschließen.

mfg
 
Zuletzt bearbeitet:

Kibo

Registriertes Mitglied
Prognosevariablen
%Planeten%.m(1, 2, 3 oder 4)v(x,y oder z)
Ist das nicht die?

Oder schlägst du vor das RK umzubauen? zum Beispiel:

einmal RK zum integrieren Beschleunigung>Geschwindigkeit und nachträglich noch einmal RK für integrieren Geschwindigkeit>Koordinaten
Oder eins von beiden eulern?
oder ganz anders?

mfg
 

Bernhard

Registriertes Mitglied
Ist das nicht die?
Ich würde die m-Variablen eher Steigungen nennen. Der Prognosepunkt und die Prognosegeschwindigkeit berechnen sich dann gemäß:

Px = x + Schrittweite / 2 * m1v
Pv = v + Schrittweite / 2 * m1a

usw., denn es gibt dann zu jedem m ein Prognosepaar.

Oder schlägst du vor das RK umzubauen?
Das Runge-Kutta-Verfahren ist fix. Es muss nur noch implementiert werden und zwar am besten fehlerfrei ;) .

Oder eins von beiden eulern?
Pfui.
 
Oben