Stimmt diese Planetenlaufbahn-Approximation?
Heyy ich versuche seit einiger Zeit mit Python eine Animation für die Laufbahn von der Erde zu erstellen, aber ich wundere mich, ob mein Ansatz stimmt. 🤔
Ich will ein Python Programm erstellen, wo ich den Exponenten des Radius im Gravitationsgesetz ändern und somit neue Umlaufbahnen zeichnen kann. Dabei stütze ich mich auf die Formel: F = G * (m1 * m2)/r^x. Ich bin mir nicht sicher, ob mein Programm korrekte Laufbahnen zeichnet, da ich nicht sonderlich gut in Physik bin. Der Grundgedanke kann ich nachvollziehen: Bei abnehmender Kraft, wird die Umlaufbahn grösser. Aber ist es normal, dass die Umlaufbahn der Erde bei einer kleinen Exponenten-Änderung von 2 auf 2.02 so viel grösser wird? Falls ja, warum ist das so?
Danke im Voraus für eure Hilfe! 😊
LG Leonie
Ein Quadrat entspricht hier einer astronomischen Einheit.
Hier ist mein Programm:
Ich verwende für meine Animation pygame.
import pygame
import sys
import math
Das sind die Konstanten, die ich verwende und die Startwerte:
G = 6.67430e-11
M_SUN = 1.989e30
M_EARTH = 5.972e24
AU = 1.496e11
TIME_STEP = 43200
x = AU
y = 0
vx = 0
vy = 29780
exp_r = float(input"Exponent:"))
Die Erstellung des Fensters:
pygame.init()
WIDTH, HEIGHT = 800, 800
SCALE = WIDTH / (20 * AU) # Maßstab für die Darstellung
screen = pygame.display.set_mode((WIDTH, HEIGHT))
clock = pygame.time.Clock()
Mit RK4 approximiere ich die Laufbahn der Erde. Hier sind die Funktionen:
def runge_kutta4(t, state, dt):
k1 = deriv(t, state)
k2 = deriv(t + 0.5 * dt, [state[i] + 0.5 * dt * k1[i] for i in range(4)])
k3 = deriv(t + 0.5 * dt, [state[i] + 0.5 * dt * k2[i] for i in range(4)])
k4 = deriv(t + dt, [state[i] + dt * k3[i] for i in range(4)])
return [state[i] + (dt / 6) * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) for i in range(4)]
def draw_grid():
grid_size = AU * SCALE
half_width = WIDTH // 2
half_height = HEIGHT // 2
for i in range(-19, 20):
x_pos = int(half_width + i * grid_size)
pygame.draw.line(screen, (50, 50, 50), (x_pos, 0), (x_pos, HEIGHT))
for i in range(-19, 20):
y_pos = int(half_height + i * grid_size)
pygame.draw.line(screen, (50, 50, 50), (0, y_pos), (WIDTH, y_pos))
Hauptschleife:
running = True
t = 0
earth_path = []
while running:
for event in pygame.event.get():
if event.type == pygame.QUIT:
running = False
state = [x, y, vx, vy]
x, y, vx, vy = runge_kutta4(t, state, TIME_STEP)
t += TIME_STEP
draw_x = int(WIDTH / 2 + x * SCALE)
draw_y = int(HEIGHT / 2 + y * SCALE)
earth_path.append((draw_x, draw_y))
screen.fill((0, 0, 0))
draw_grid()
pygame.draw.circle(screen, (255, 255, 0), (WIDTH // 2, HEIGHT // 2), 20)
pygame.draw.circle(screen, (0, 0, 255), (draw_x, draw_y), 5)
for px, py in earth_path:
pygame.draw.circle(screen, (0, 255, 0), (px, py), 1)
pygame.display.flip()
clock.tick(60)
pygame.quit()
sys.exit()
Hier ist noch die deriv-Funktion:
def deriv(t, state):
x, y, vx, vy = state
x, y, vx, vy = Decimal(x), Decimal(y), Decimal(vx), Decimal(vy)
r = (x**2 + y**2).sqrt()
if r == 0:
return [vx, vy, Decimal("0"), Decimal("0")]
r_modified = r**exp_r
F = -G * M_SUN / r_modified
ax = F * x / r
ay = F * y / r
return [vx, vy, ax, ay]
3 Antworten
Ich habe mir Deine Runge-Kutta-Funktion nicht näher angeschaut (und schon garnicht nachgerechnet).
Mir ist jedoch aufgefallen, das diese mit Floats nahe den Min/Max-Values operiert.
1.175494351E-38 und 3.402823466E+38
Float benutzt 32 Bit und "rundetc um so "aggressiver", je näher Resultate von Operationen an den Limits liegen. Dies kann zu massiven Rechenfehler führen, da jeweils nur eine stark begrenzte Anzahl an Bits für Exponent und Mantisse zur Verfügung stehen.
Versuche es statt dessen mal mit Double oder Decimal oder sogar BigInt (sollte auch Python haben).
Double rechnet mit 64Bit , Decimal mit 128 und Bigint sogar mit ca. max. 20Mrd Bit (extrem langsam auf Konsumerrechnern).
Für BigInt-Berechnungen musst Du gegebenenfalls kommazahlen manuell in (big)Integer scalieren.
Ich habe zu dieser frühen Stunde leider keine Zeit auf Deine Berechnungen näher einzugehen. Versuche es erstmal mit den übrigen (genaueren) Datentypen zu arbeiten. Sollte dies nicht klappen werde ich mich mal über Deine Funktion beugen.
Python kennt nur einen eingebauten Gleitkommatyp, und der wird auf allen modernen Plattformen in 64 Bit untergebracht, entspricht also "double".
Vielen Dank für deine Hilfe! 😊 Ich nun die Dezimal-Funktion in meinem Programm eingefügt, welche deutlich mehr Nachkommastellen ermöglicht. Nun werden die Ellipsen zwar etwas schöner gezeichnet, aber deren Form und Grösse oder besser gesagt deren Exzentrizität bleibt gleich...
Aber ist es normal, dass die Umlaufbahn der Erde bei einer kleinen Exponenten-Änderung von 2 auf 2.02 so viel grösser wird? Falls ja, warum ist das so?
Ja, das sollte so schon richtig sein. Der Exponent bestimmt, wie stark die Anziehungskraft mit der Entfernung abnimmt. Bei Werten über 2 (Normalwert) nimmt diese Kraft stärker ab, also wird die Erde weniger zurückgezogen und eine größere, ovale Umlaufbahn entsteht. Das macht sich schon bei kleineren Zuwächsen deutlich bemerkbar. Deshalb wird diese Kraft auch mit einer Potenzfunktion abgebildet. Der Wert steigt (oder sinkt) exponentiell.
Probiere für deine Simulation doch noch weitere Werte. Womöglich geht es bei einem Exponent wie 2.5 bereits in eine Hyperbelbahn.
Hier ist mein Programm:
1) Deine deriv-Funktion hast du gar nicht geteilt, dabei sollte die doch erst aufzeigen, ob du exp_r richtig verwendest.
2) Hier fehlt die Klammer vor dem ersten Anführungszeichen:
exp_r = float(input"Exponent:"))
3) Auch wenn du deine mathematischen Formeln als Vorlage hast, würde ich dir empfehlen, Bezeichner wie AU, dt, k1, t, u.ä. gegen aussagekräftigere Namen zu tauschen.
4) Für mehr (wissenschaftliche) Genauigkeit würde ich nicht PyGame für die Visualisierung nutzen, sondern eher Matplotlib oder vielleicht Manim.
Die Konversion der Werte in decimal kommt zu spät, wenn du sie in der deriv-Funktion behältst. Die Werte für x, y, usw. werden ja viel früher, in runge_kutta4 berechnet, wo Python noch seinen float-Typ benutzen wird. Bei den Berechnungen müsstest du darauf achten, dass alle Fließkommazahlterme explizit in decimal angelegt/konvertiert werden, denn sonst kann es zu einem Typfehler kommen (decimal und float können nicht direkt miteinander verrechnet werden).
Meiner Meinung reicht es aber aus, für diese Simulation nur mit float zu arbeiten. Wenn es bei deiner Rechnung zu (physikalischen) Ungenauigkeiten kommt, dann sind da andere Quellen bedeutender (z.B. das Annäherungsverfahren via Runge-Kutta oder die Präzision der Konstanten), wobei ich auch da nicht andeuten möchte, dass du dir eine andere Lösung suchen solltest - dein Lösungsweg ist vollkommen legitim.
Du könntest zur Prüfung noch die jeweilige Geschwindigkeit der Erde berechnen und mit ihrer Fluchtgeschwindigkeit vergleichen.
escape_velocity = ((2 * G * M_SUN) / r).sqrt()
Sobald diese Geschwindigkeit erreicht wird, fällt die Erde aus ihrer Bahn.
Danke für deinen Tipp! 😃 Ich habe jetzt zu Beginn der Funktion Runge Kutta folgende Zeile hinzugefügt: state = [Decimal(s) for s in state] . Bei einer Geschwindigkeit von 43 km/s bei einem Exponenten r von 2, fliegt die Erde tatsächlich aus dem Sonnensystem heraus. 👍
Mich überrascht das nicht: Beim Exponenten 2.02 hast Du nur noch 60% Anziehungskraft, bei 2.03 nur noch 46%. Kein Wunder, dass sich die Erde da aus dem Staub macht.
Gib ihr wenigstens eine Chance auf eine stabile Umlaufbahn, indem Du sie langsamer machst, näher an die Sonne rückst oder G erhöhst.
Vielen Dank für deinen Kommentar! 👍 Es ist erstaunlich wie schnell die Anziehungskraft abnimmt, das hätte ich zu Beginn des Programmierens nicht gedacht. 😅
So schnell wird das nicht weniger. Zu Hause dürftest Du keinen Unterschied merken. Aber auf 150 Milliarden Metern läppert sich doch etwas zusammen.
Danke für deine Hilfe! Ich habe nun die deriv-Funktion bei meiner Fragestellung hinzugefügt und werde auch nach deinem Tipp ein Matplotlib-Diagramm erstellen. 😊 Mich wundert es immer noch, wie schnell sich die Laufbahn ändert. Ab einem Exponent-Wert r von 2.03 wird gar keine Ellipse mehr gezeichnet: Sie wird nur angefangen und anschliessend wird die Erde aus dem Sonnensystem geschleudert. Bei 2.5 wird sogar nur noch eine gerade, senkrechte Linie nach unten gezeichnet. 😆