Gravitationssimulation - Umlaufbahn keine Ellipse

Wissenschaft Physik
1

Hallo,

ich möchte mit Visual Basic 6 eine Gravitationssimulation schreiben. Im Modell gibt es zunächst in einem zweidimensionalen Raum einen Planeten und ein Objekt,das von dem Planet angezogen wird, und dessen Orbit ich mit der Simulation darstellen will.

Dazu habe ich folgende Rechnung angestellt:

Variablen (Vektoren sind großgeschrieben, Zahlen klein):

O, P: Ortsvektoren des Objektes und des Planeten
R: Richtungsvektor, der angibt in welche Richtung die Gravitation wirkt (Einheitsvektor)
A: Momentane Gravitaionsbeschleunigung

k: Konstante
f: Gravitaionskraft
mp: Masse des Planeten
mo: Masse des Objektes

Zuerst wird der Richtungsvektor R errechnet und zum Einheitsvektor gemacht, die länge wird später mit dem Gravitationsgesetz errechnet.

R = (P-O) / Betrag(P-O)

Danach nach dem Gravitaionsgesetz f errechnet

f = k * mp * mo / (Betrag(P-O)²)

Nun wird die Gravitationsbeschleunigung A nach F = m*a (A = F / m) mit Richtung und entsprechender Länge errechnet

A = (R * f) / m

danach wird A nur noch zur momentanen Geschwindigkeit addiert und der Ortsvektor O mit der Geschwindigkeit verrechnet.

Diese Berechnungen erfolgen immer im gleichen Zeitabstand (ca. 10 ms). Merkwürdigerweise ist die entstehende Umlaufbahn keine Ellipse wie sie nach Kepler eigentlich sein sollte, sondern eine sich ständig ändernde Umlaufbahn; Egal welche Anfangswerte ich für die Geschwindigkeit wähle (und das Objekt natürlich nicht kollidiert oder aus dem Gravitationsfeld kommt).
Ich dachte zuerst dass es an der hier vernächlässigten Gegenwirkung (Gravitationswirkung von dem Objekt auf den Planeten) liegt, als ich diese hinzugefügt habe hat das aber nichts geändert.
Rundungsfehler kann ich eigentlich auch ausschließen da ich den Datentyp Double verwende, der sehr genau ist.

Wo liegt der Fehler in meiner Rechnung?

mfg
Raber

2

Möglichkeiten
Moin,

ohne Zahlenwerte zu haben, ist es schwierig zu sagen, wo die Fehler stecken könnten. Zwei Stellen, wo sie möglicherweise sind, sind mir aufgefallen:

Nun wird die Gravitationsbeschleunigung A nach F = m*a (A = F
/ m) mit Richtung und entsprechender Länge errechnet

A = (R * f) / m

Hast Du aufgepaßt, daß Du hier die richtige Masse einsetzt? Also mo?

danach wird A nur noch zur momentanen Geschwindigkeit addiert

Vorher muß A noch mit der Zeit multipliziert werden gemäß v=a*t. Hast Du das gemacht?

Mehr fällt mir ohne konkrete Zahlen auch nicht ein.

Gruß Kubi

3

Hallo Raber,

Wo liegt der Fehler in meiner Rechnung?

Deine Überlegungen sind grundsätzlich richtig; ich habe darin jedenfalls nichts grob Fehlerhaftes finden können. Aber vielleicht hat sich irgendein subtiler Bug in den Code eingeschlichen, den Du nicht erkennst – so einer von denen, die an sich klein sind, aber eine verheerende Wirkung haben (z. B. Pluszeichen statt Minuszeichen o. ä.). Vorschlag: Poste doch einfach mal die entsprechende VB-Routine hier.

Mit freundlichem Gruß
Martin

4

Hallo,

wieso nimmst du nicht gleich die Ellipsengleichung? Dann hast du direkt die exakte Bahn. Diese Interpoliererei ist doch immer fehleranfällig.

Gruß
Oliver

5

wieso nimmst du nicht gleich die Ellipsengleichung?

Hallo Oliver,

verwendet man nur die Bahngleichung, dann hat man nicht die Dynamik, man weiß also nicht, zu welchem Zeitpunkt der Planet an welcher Stelle ist. Man weiß nur, dass der Planet irgendwann mal jede Stelle der Bahnkurve durchläuft.

Viele Grüße
Stefan

6

Nun wird die Gravitationsbeschleunigung A nach F = m*a (A = F
/ m) mit Richtung und entsprechender Länge errechnet

A = (R * f) / m

danach wird A nur noch zur momentanen Geschwindigkeit addiert
und der Ortsvektor O mit der Geschwindigkeit verrechnet.

Diese Berechnungen erfolgen immer im gleichen Zeitabstand (ca.
10 ms).

Hallo Fredde,

zusätzlich zu den bereits genannten Fehlermöglichkeiten fällt mir noch folgendes ein:
Genauer wird Deine Rechnung natürlich, wenn Du die 10ms kleiner werden lässt, aber diese Rechnung bleibt eine Näherung, es sei denn, diese kurze Zeitspanne wird unendlich klein. Dann musst Du allerdings unendlich lange auf Dein Ergebnis warten.

Das von Dir benutzte Verfahren nennt sich Vorwärts-Euler-Verfahren. Es ist das einfachste Verfahren zur numerischen Lösung von Differentialgleichungssystemen. Eines mit besserer Konvergenz ist zum Beispiel das Runge-Kutta-Verfahren.
So, das waren jetzt einige Begriffe, mit denen man Google füttern kann. Ich wünsche Dir noch viel Spaß bei weiteren Berechnungen. Du dringst in einen äußerst interessanten Themenbereich ein.

Viele Grüße
Stefan

7

Hallo,

ich möchte mit Visual Basic 6 eine Gravitationssimulation
schreiben. …
danach wird A nur noch zur momentanen Geschwindigkeit addiert
und der Ortsvektor O mit der Geschwindigkeit verrechnet.

Hallo,

du rechnest also die Beschleunigung und die Geschwindigkeit zum Zeitpunkt t, t+10ms, t+20ms usw. Das ergibt natürlich einen systematischen Fehler, da du die Geschwindigkeit als konstant für 10 ms annimmst, was nicht stimmt.

Die (korrekte) Bahn besteht aus 2 Hälften, in der einen wird der Planet ständig beschleunigt, in der anderen ständig abgebremst. Deine Rechenmethode liefert also für die erstere immer zu kleine Geschwindigkeitswerte, für die andere Hälfte immer zu grosse.

Wenn ich das ganz grob nach Gefühl interpretiere, müsste deine Bahn von der wahren Bahn auf der beschleunigten Hälfte nach innen abweichen und auf der gebremsten Hälfte nach aussen, es ist kein Wunder, wenn sich keine geschlossene Bahn mehr ergibt.

Das ganze sieht dann eher aus wie die relativistische Periheldrehung des Merkur, aber die hast du ganz bestimmt nicht berechnet.

Gruss Reinhard

8

Hallo,

deine Antwort erscheint mir ziemlich plausibel, tatsächlich weichen die Werte auch wie beschrieben ab.
Der Fehler liegt warscheinlich dann darin dass in den Formeln kein t auftaucht. Aber wie bekomme ich das dort hinein?

Wie gesagt wird alle 10ms die Richtung und die Kraft neu berechnet und diese dann zur momentanen Geschwindigkeit hinzuaddiert. Mithilfe der Gravitationsformel kann ich ja eben nur die momentan wirkende Kraft errechnen.
Gibt es denn eine Formel mit der man ähnlich wie beim freien Fall s(t) anhand von Masse, Koordinaten und Anfangsgeschwindigkeiten berechnen kann?

mfg
Raber

[Bei dieser Antwort wurde das Vollzitat nachträglich automatisiert entfernt]

9

Hi,

hier ist der gewünschte Quelltext mit den entsprechenden Werten

Wenn ich allerdings die Timerzeit verkürze ändert das nichts an der Rechnung, weil ja meine Berechnungen Zeitunabhängig sind, schließlich taucht in den Formeln ja nirgends t auf, oder?
Es wird nach jedem Zeitintervall anhand der Koordinaten die momentane Kraft und Beschleunigung ausgerechnet, die auf den Körper wirkt, und dann diese Beschleunigung zur Bewegung hinzuaddiert, welche Wiederum die Position des Körpers geringfügig ändert. Im nächsten Zyklus wird dasselbe mit den neuen Koordinaten gemacht. Die Genauigkeit der Rechnung ist daher unabhängig von der Länge des Intervalls.
Daher wird A ja auch nicht mit t multipliziert bevor es zur Geschwindigkeit addiert wird, weil sich der Richtungsvektor des Objekts ja ständig ändert. Somit ist t immer die Länge des Intervalls (10ms) und dadurch eine Konstante. Gibt es da etwa eine Formel mit der man den Ortsvektor vom Objekt nach einer bestimmten Zeit errechnen kann (ich kenne da keine)?

Wenn jemand den Compiler für VB6 hat, kann ich dem auch gerne das ganze Projekt zuschicken, dann wird das noch anschaulicher.

mfg
Raber

Public Type MathVector
X1 As Double
X2 As Double
X3 As Double
End Type

Private Type Planet
Ort As MathVector
Beschleunigung As MathVector
Masse As Double
End Type

Private Sub Gravitation(Angezogener As Planet, Anzieher As Planet)

Dim RP As MathVector
Dim Richtung As MathVector
Dim Gravi As MathVector
Dim u As Double

Const f = 10000 'Gravitationskonstante

RP = VectorMinus(Anzieher.Ort, Angezogener.Ort)

'Einheitsvektor herstellen/Richtung berechnen
Richtung = SMulti(1 / VectorBetrag(RP), RP)

'Vektorlänge errechnen
u = -(f * Angezogener.Masse * Anzieher.Masse) / (VectorBetrag(RP) ^ 2)

'Gravitation gesamt berechnen
Gravi = SMulti(u, Richtung) 'Gravitationsfaktor
Gravi = SMulti(1 / Angezogener.Masse, Gravi) 'Massefaktor

Angezogener.Beschleunigung = VectorPlus(Angezogener.Beschleunigung, Gravi)

Angezogener.Ort = VectorPlus(Angezogener.Ort, Angezogener.Beschleunigung)

End Sub

'Berechnet S-Multiplikation
Function SMulti(ByVal Skalar As Double, Vektor As MathVector) As MathVector
SMulti.X1 = Vektor.X1 * Skalar
SMulti.X2 = Vektor.X2 * Skalar
SMulti.X3 = Vektor.X3 * Skalar
End Function

Function VectorMinus(Vektor1 As MathVector, Vektor2 As MathVector) As MathVector
VectorMinus.X1 = Vektor2.X1 - Vektor1.X1
VectorMinus.X2 = Vektor2.X2 - Vektor1.X2
VectorMinus.X3 = Vektor2.X3 - Vektor1.X3
End Function

Function VectorPlus(Vektor1 As MathVector, Vektor2 As MathVector) As MathVector
VectorPlus.X1 = Vektor2.X1 + Vektor1.X1
VectorPlus.X2 = Vektor2.X2 + Vektor1.X2
VectorPlus.X3 = Vektor2.X3 + Vektor1.X3
End Function

Function VectorBetrag(Vektor As MathVector) As Long
VectorBetrag = Sqr(Vektor.X1 ^ 2 + Vektor.X2 ^ 2 + Vektor.X3 ^ 2)
End Function

Zu Testzwecken habe ich folgende Anfangswerte eingesetzt:

Rakete.Ort.X1 = 4050
Rakete.Ort.X2 = 1000
Rakete.Beschleunigung.X1 = 200
Rakete.Masse = 5

Planet1.Ort.X1 = 5000
Planet1.Ort.X2 = 4000
Planet1.Masse = 20000

10

Auch hi,

Wenn ich allerdings die Timerzeit verkürze ändert das nichts
an der Rechnung, weil ja meine Berechnungen Zeitunabhängig
sind, schließlich taucht in den Formeln ja nirgends t auf,
oder?

so ist es!

Es wird nach jedem Zeitintervall anhand der Koordinaten die
momentane Kraft und Beschleunigung ausgerechnet, die auf den
Körper wirkt, und dann diese Beschleunigung zur Bewegung
hinzuaddiert, welche Wiederum die Position des Körpers
geringfügig ändert. Im nächsten Zyklus wird dasselbe mit den
neuen Koordinaten gemacht. Die Genauigkeit der Rechnung ist
daher unabhängig von der Länge des Intervalls.

ja! Ich bin mir nicht hundertprozentig sicher, aber eventuell ist das sogar die Ursache für die fehlerhafte Ausgabe Deines Programms.

Daher wird A ja auch nicht mit t multipliziert bevor es zur
Geschwindigkeit addiert wird, weil sich der Richtungsvektor
des Objekts ja ständig ändert. Somit ist t immer die Länge des
Intervalls (10ms) und dadurch eine Konstante. Gibt es da etwa
eine Formel mit der man den Ortsvektor vom Objekt nach einer
bestimmten Zeit errechnen kann (ich kenne da keine)?

Du mußt unterscheiden zwischen t und dt. t ist die Zeit, die von Null ab hochläuft bis z. B. 5 Minuten. Dabei tickt die Uhr in dt-Schritten. dt kann z. B. 0.1 s lang sein; dann wäre die Simulation nach 3000 Zeitschritten fertiggerechnet (dann sind die 5 min erreicht).

Die Formeln lauten:

[fertig berechnet an dieser Stelle sind ax und ay]

// "Update" der Geschwindigkeit v\>
vx := vx + ax\*dt;
vy := vy + ay\*dt;

// "Update" der Position r\>
rx := rx + vx\*dt
ry := ry + vy\*dt

// Uhr um einen Zeitschritt weiterticken lassen
t := t + dt

Du solltest das Deinem Code auf jeden Fall entsprechend abändern, damit Du auch die Zeit in Deiner Simulation „drinhast“. Du kannst dann auch mal dt variieren, und herausfinden, wie die Bahnen zunehmend „falscher“ werden, wenn Du dt zu groß wählst.

Und noch ein Tip: Verwende für die Gravitationskonstante doch den echten Wert (6.672E-11). Dann kannst Du für die Massen auch echte Werte einsetzen (z. B. Masse der Erde, Masse eines Satelliten), und für alle anderen Simulationsparameter (Anfangsbedingungen, Zeitschritt dt) ebenfalls. Dann kannst Du Dir angucken (grafische Ausgabe vorausgesetzt), wie ein z. B. ein Projektil bei einer Abschußgeschwindigkeit von 11.2 km/s (Fluchtgeschwindigkeit der Erde) tatsächlich auf Nimmerwiedersehen im All verschwindet (gerade so) und ähnliche „Scherze“.

Gruß
Martin

11

Hallo,

ich habe jetzt den Zeitfaktor eingebaut, zur Beschleunigung und Geschwindigkeit wird per S-Multiplikation noch der Faktor dt hinzugerechnet.
Leider hat es die fehlerhafte Flugbahn nicht geändert, jedenfalls nicht auffallend.
Ich habe die Berechnung in eine Endlosschleife gesetzt und kann dadurch noch gerade so mit einem dt von 0,001 rechnen, ansonsten wäre mein Rechner überfordert. Die Umlaufbahn ist leider exakt dieselbe.

Warscheinlich müsste ich den Zeitabstand unendlich klein machen, also den Grenzwert gegen null bilden. Aber wie kann ich das hier machen, schließlich ist es ja nicht eine einfache Funktion die abgeleitet werden kann (sehr viele Parameter!)?
Wie das Runge-Kutta-Verfahren funktioniert kann ich leider nicht verstehen…

mfg
Raber

[Bei dieser Antwort wurde das Vollzitat nachträglich automatisiert entfernt]

12

Hallo,

ich habe jetzt den Zeitfaktor eingebaut, zur Beschleunigung
und Geschwindigkeit wird per S-Multiplikation noch der Faktor
dt hinzugerechnet.
Leider hat es die fehlerhafte Flugbahn nicht geändert,
jedenfalls nicht auffallend.

hmm… poste den neuen Code!

Runge-Kutta ist toll, brauchst Du aber für dieses Projekt mit Sicherheit nicht. Wenn alles korrekt ist und dt weder viel zu groß noch viel zu klein gewählt ist, stehen die Bahnkurven wie ne Eins.

Gruß
Martin

13

Gibt es denn eine Formel mit der man ähnlich wie beim freien
Fall s(t) anhand von Masse, Koordinaten und
Anfangsgeschwindigkeiten berechnen kann?

Eine Möglichkeit gäbe es da noch, die dir garantiert eine Ellipse liefert.

Und zwar gehst du aus von der Gleichung für die Ellipse:

r = r(phi)

die den Radius in Abhängikeit des Winkels angibt.

Dann gibt es noch eine Integralgleichung von der Art

t® = Integral(r0 bis r)dr’ f(r’)

Sie gibt die Zeit in Abhängigkeit des Radiuses an.

Du könntest jetzt folgendes tun:
Wenn du in obiger Integralgleichung für r0 den Wert der kleinen Halbachse a setzt, und dann r zwischen a und der großen Halbachse b variierst, hast du:

t®, für a

14

Hallo,
bei einer Simulation muß man immer betrachten, wie groß denn der Fehler des ganzen ist. Vielleicht sind Deine 10ms zu groß, vielleicht auch zu klein. Um das abzuschätzen, berechnet man bei der Simulation erstmal mit einer Schrittgröße, dann mit einer kleineren (z.B. Faktor 10 kleiner). Wenn die Differenz zwischen den beiden Versionen klein genug ist, bleibt man bei der größeren Schrittweite. Wenn die Differenz zu groß ist, wiederholt man die Rechnung wiededrum mit nochmals kleinerer Schrittweite und so weiter.
Zum Zweiten fällt mir noch ein: mit welchen Startbedingungen hast Du denn angefangen? Die Planeten wurden ja auch nicht einfach stehend im All platziert, sondern sind an ‚passender‘ Stelle mit bereits ‚passendem‘ Energie- und Bewegungsvektor entstanden.
Gruß
Axel

15

Hallo,

ich habe gestern bei den Tests vergessen, dass ich ja ein Doppelsternsystem habe, dadurch war der Orbit keine Ellipse (ich dachte eigentlich das wäre auch dort so). Ich war eben schon etwas müde … Tut mir leid dass ich es nicht früher bemerkt habe.

Jedenfalls ist bei den Berechnungen mit einem Planeten für ein kleines dt tatsächlich eine (fast-)Ellipse rausgekommen. Offenbar war die Grobheit also der Fehler. Bleibt allerdings immer noch das Problem, dass ein kleinerer Berechnungsfehler immer übrig bleibt. Auch wenn er evtl. zu vernächlässigen wäre möchte ich doch einen exakten Simulator programmieren; schließlich bin ich ja Perfektionist.
Wie ginge das denn mit dem Runge-Kutta-Verfahren?

mfg
Raber

[Bei dieser Antwort wurde das Vollzitat nachträglich automatisiert entfernt]

16

Hallo Oliver,

Gibt es denn eine Formel mit der man ähnlich wie beim freien
Fall s(t) anhand von Masse, Koordinaten und
Anfangsgeschwindigkeiten berechnen kann?

Eine Möglichkeit gäbe es da noch, die dir garantiert eine
Ellipse liefert.

Und zwar gehst du aus von der Gleichung für die Ellipse:

r = r(phi)

die den Radius in Abhängikeit des Winkels angibt.

Legst Du den Ursprung des Polarkoordinatensystems in den Mittelpunkt oder in einen Brennpunkt der Ellipse?

Dann gibt es noch eine Integralgleichung von der Art

t® = Integral(r0 bis r)dr’ f(r’)

Sie gibt die Zeit in Abhängigkeit des Radiuses an.

Soweit ich mich erinnere, ist dieses Integral nicht geschlossen zu lösen.

Du könntest jetzt folgendes tun:
Wenn du in obiger Integralgleichung für r0 den Wert der
kleinen Halbachse a setzt, und dann r zwischen a und der
großen Halbachse b variierst, hast du:

t®, für a

17

möchte ich doch einen exakten Simulator programmieren;
schließlich bin ich ja Perfektionist.
Wie ginge das denn mit dem Runge-Kutta-Verfahren?

mfg
Raber

Hallo Raber,

nun doch eine dringende Buchempfehlung:
Andreas Guthmann: „Einführung in de Himmelsmechanik und Ephemerdenrechnung, Theorie, Algorithmen, Numerik“ Spektrum Akademischer Verlag, ISBN3827405742

Sowohl Zweikörperprobleme, die zu den klassischen keplerschen Ellipsen führen, werden behandelt, als auch Mehrkörperprobleme, die dann nur noch mit numerischer Integration zu machen sind. Auch das Runge-Kutta-Verfahren wird erklärt. Und noch viel mehr.

Viele Grüße
Stefan

18

Hallo,

Legst Du den Ursprung des Polarkoordinatensystems in den
Mittelpunkt oder in einen Brennpunkt der Ellipse?

Der Mittelpunkt ist der Massenschwerpunkt des Systems.

Soweit ich mich erinnere, ist dieses Integral nicht
geschlossen zu lösen.

Da erinnerst du dich richtig. Das muss man dann numerisch berechnen. Dürfte aber auch kein Problem sein.

Kennst Du die Keplerschen Gesetze? Da finde ich allenfalls die
eine Hälfte der Bahn symmetrisch zur anderen Hälfte. Wie
kommst Du auf Viertel?

Da hast du natürlich recht. Man muss das zweite Viertel analog ausrechnen. Die andere Hälfte der Bahn ergibt sich dann aber.

Dem stimme ich zu. Allerdings hat diese Ellipsenbahn nichts
mit Planetenbewegungen zu tun.

Natürlich hat sie was mit der Planetenbahn zu tun: Die eingangs verwendete Ellipsengleichung r(phi) ist ja gerade die Lösung des Keplerproblems. Man beschafft sich nur noch die Dynamik entlang der Bahn. Der Vorteil ist, dass Näherungen bei jener Berechnung nichts mehr an der Bahn selbst ändern.

Gruß
Oliver