Newtons Gravitation mit GNUPlot und Python
Nach dem ich im Allgemeine Fragen Forum schon ein bischen gepostet hab, hier mein Programm, das die Bewgung beliebig vieler Massepunkte über GNUPlot visualisieren soll. Ich hoffe die Formeln stimmen und sind richtig umgesetzt, die Bezeichungen sind erklärend genug - für Physiker sind sie es... immernoch, da bin ich fest von überzeugt - und ihr habt GNUPlot!
Das Programm selbst ist eine Übersetzung eines alten PASCAL Programmes von mir und wartet deswegen nicht gerade mit pythonscher Eleganz auf. Sicherlich könnte man die Daten der Planeten alle in einer Liste speichern, aber ich denke, dass würde alles sehr unübersichtlich machen. Ein bischen was ist schon verbessert: An statt, wie in der orginal Version die Daten direkt auf den Bildschirn auszugeben, speichert das Skript die Bewegung jedes Massepunktes in einer Datei und übergibt diese an GNUPlot, womit ich mich um die Programmierung eines GUI drücken kann.
Es existieren zwei Versionen, in der ersten muss man die Daten im Quelltext selbst ändern, in der zweiten läuft alles über raw_input, auchnicht elegant, aber zu mehr sehe ich mich nicht im Stande.
Es folgt Version 1:
[Code ausgelagert]
Edit (Leonidas): Code ausgelagert.
Das Programm selbst ist eine Übersetzung eines alten PASCAL Programmes von mir und wartet deswegen nicht gerade mit pythonscher Eleganz auf. Sicherlich könnte man die Daten der Planeten alle in einer Liste speichern, aber ich denke, dass würde alles sehr unübersichtlich machen. Ein bischen was ist schon verbessert: An statt, wie in der orginal Version die Daten direkt auf den Bildschirn auszugeben, speichert das Skript die Bewegung jedes Massepunktes in einer Datei und übergibt diese an GNUPlot, womit ich mich um die Programmierung eines GUI drücken kann.
Es existieren zwei Versionen, in der ersten muss man die Daten im Quelltext selbst ändern, in der zweiten läuft alles über raw_input, auchnicht elegant, aber zu mehr sehe ich mich nicht im Stande.
Es folgt Version 1:
[Code ausgelagert]
Edit (Leonidas): Code ausgelagert.
Und hier, der Übersichthalber in einem zweiten Post, die zweite Variante:
[Code ausgelagert]
Edit (Leonidas): Code ausgelagert.
[Code ausgelagert]
Edit (Leonidas): Code ausgelagert.
- birkenfeld
- Python-Forum Veteran
- Beiträge: 1603
- Registriert: Montag 20. März 2006, 15:29
- Wohnort: Die aufstrebende Universitätsstadt bei München
Tu bitte Lesern deines Codes den Gefallen und achte auf die Zeilenlänge. Gerade Kommentare kann man sehr schön in die Zeile vor dem kommentierten Code setzen... wenn man schon bei maximiertem Fenster eine horizontale Scrollbar braucht, ist das unschön.
Außerdem wäre die Option "-persist" für Gnuplot angebracht, damit das Plotfenster nicht gleich wieder verschwindet.
Außerdem wäre die Option "-persist" für Gnuplot angebracht, damit das Plotfenster nicht gleich wieder verschwindet.
Die ganzen
sind sehr unpythonisch. Besser wäre z.B. für
oder
Die while True-Schleife könnte man besser als
formulieren.
Die Kraft zwischen 2 Planeten A und B wird doppelt berücksichtigt, einmal von A nach B, dann von B nach A. Wenn man zur Korrektur die Kraft halbiert, ändert sich durchaus auch die Form der Kurve.
Das Programm ließe sich sicher ideal in Numpy umsetzen, was wohl schneller wäre und auch die Erweiterung auf die dritte Dimension einfach machen würde.
MfG
HWK
Code: Alles auswählen
for i in range(anzplan):
Code: Alles auswählen
for i in range(anzplan):
f[i][0] = 0
f[i][1] = 0
Code: Alles auswählen
for f_ in f:
f_[0] = 0
f_[1] = 0
Code: Alles auswählen
for f_ in f:
f_[:] = [0, 0]
Code: Alles auswählen
for count in range(interationen):
Die Kraft zwischen 2 Planeten A und B wird doppelt berücksichtigt, einmal von A nach B, dann von B nach A. Wenn man zur Korrektur die Kraft halbiert, ändert sich durchaus auch die Form der Kurve.
Das Programm ließe sich sicher ideal in Numpy umsetzen, was wohl schneller wäre und auch die Erweiterung auf die dritte Dimension einfach machen würde.
MfG
HWK
Ernsthaft? Ich hab mit den Kräften so wie so noch Probleme, in der neusten Version werden alle 50 Durchläufe Kraftvektoren eingezeichnet, aber die stimmen hinten und vorne nicht...HWK hat geschrieben:Die Kraft zwischen 2 Planeten A und B wird doppelt berücksichtigt, einmal von A nach B, dann von B nach A. Wenn man zur Korrektur die Kraft halbiert, ändert sich durchaus auch die Form der Kurve.
http://paste.pocoo.org/show/3986/
- birkenfeld
- Python-Forum Veteran
- Beiträge: 1603
- Registriert: Montag 20. März 2006, 15:29
- Wohnort: Die aufstrebende Universitätsstadt bei München
Das ist vollkommen in Ordnung. Beide Planeten wechselwirken und spüren beide eine Kraft (siehe Newtons zweites Gesetz).HWK hat geschrieben: Die Kraft zwischen 2 Planeten A und B wird doppelt berücksichtigt, einmal von A nach B, dann von B nach A.
- birkenfeld
- Python-Forum Veteran
- Beiträge: 1603
- Registriert: Montag 20. März 2006, 15:29
- Wohnort: Die aufstrebende Universitätsstadt bei München
Möglicherweise liegt das daran, dassAgroschim hat geschrieben: Ernsthaft? Ich hab mit den Kräften so wie so noch Probleme, in der neusten Version werden alle 50 Durchläufe Kraftvektoren eingezeichnet, aber die stimmen hinten und vorne nicht...
Code: Alles auswählen
r[i]**1/2
Für die Berechnung des Kraftvektors würde ich übrigens direkt die vektorielle Formel hernehmen:
Code: Alles auswählen
gamma*m1*m2
F = -------------- ( x1-x2, y1-y2 )
|r_12|^3
Probleme gelöst, hier ist der Code:
http://paste.pocoo.org/show/3994/
Und hier ein Beispiel:
http://freenet-homepage.de/agroschim/planeten.png
http://paste.pocoo.org/show/3994/
Und hier ein Beispiel:
http://freenet-homepage.de/agroschim/planeten.png
Stimmt. Da die Kräfte für beide Planeten aber bis auf das Vorzeichen identisch sind, kann man dann die Anzahl der Berechnungen halbieren:birkenfeld hat geschrieben:Das ist vollkommen in Ordnung. Beide Planeten wechselwirken und spüren beide eine Kraft (siehe Newtons zweites Gesetz).
Code: Alles auswählen
for ii in range(anzplan):
for i in range(ii + 1, anzplan):
ds0 = s[ii][0] - s[i][0]
ds1 = s[ii][1] - s[i][1]
r = ds0 * ds0 + ds1 * ds1
if r:
fg = g * m[i] * m[ii] / r ** (3 / 2)
temp = ds0 * fg
f[i][0] += temp
f[ii][0] -= temp
temp = ds1 * fg
f[i][1] += temp
f[ii][1] -= temp
Für r und ds sind übrigens keine Listen notwendig.
MfG
HWK
Zuletzt geändert von HWK am Montag 17. September 2007, 20:44, insgesamt 1-mal geändert.
Code: Alles auswählen
a b c d
a x ab ac ad
b -ab x bc bd
c -ac -bc x cd
d -ad -bd -cd x
Code: Alles auswählen
fg / r[i] ** (1 / 2)
Code: Alles auswählen
fg = g * m[i] * m[ii] / r ** (3 / 2)
Code: Alles auswählen
f[i][0] += temp
Code: Alles auswählen
f[i][0] = f[i][0] + temp
Code: Alles auswählen
ds0 * ds0
Code: Alles auswählen
ds0 ** 2
Eine etwas korrektere und kürzere Variante für
Code: Alles auswählen
for i in range(anzplan):
fname = 'planet_' + str(i)
puffer = str(plot[i])
puffer = puffer.replace('[','')
puffer = puffer.replace('], ','\n')
puffer = puffer.replace(',','')
f = file(fname, 'w')
f.write (puffer)
f.close()
Code: Alles auswählen
for i, plot_ in enumerate(plot):
f = file('planet_%i' % i, 'w')
f.write ('\n'.join(('%s %s' % (x, y) for x, y in plot_)))
f.close()
MfG
HWK
- birkenfeld
- Python-Forum Veteran
- Beiträge: 1603
- Registriert: Montag 20. März 2006, 15:29
- Wohnort: Die aufstrebende Universitätsstadt bei München
Und wenn man sich schnell eine Vektor- und eine Planetenklasse bastelt, dann sehen die Berechnungen auch gleich viel schöner aus:
http://paste.pocoo.org/show/4006/ (ungetestet)
http://paste.pocoo.org/show/4006/ (ungetestet)
Zuletzt geändert von birkenfeld am Dienstag 18. September 2007, 08:03, insgesamt 1-mal geändert.
Ja, da erkennt man den Profi. Das sieht toll aus. So richtig läuft es aber nicht. Z.B. ist __div__ für Vec2 nicht definiert und count heißt jetzt iteration. Aber auch nach diesen Korrekturen sieht die Grafik noch nicht richtig aus. Das muss man sich wohl noch mal im Detail anschauen.
In Agroschims Script funktioniert übrigens die Vektorausgabe nicht richtig. Ich hab's so geändert:MfG
HWK
In Agroschims Script funktioniert übrigens die Vektorausgabe nicht richtig. Ich hab's so geändert:
Code: Alles auswählen
#Punkte plotten
for i in range(anzplan):
kommando += " 'planet_%i' w l" % i
#Vektoren plotten, bei Bedarf auskommentieren oder hinzufügen
kommando += ", 'vektoren_%i' with vectors head filled lt 2" % i
if i != anzplan - 1:
kommando += ", "
HWK
@birkenfeld: Das Problem lag an einer falschen Formel von Agroschim. Bei ihm und bei mir klappte es zufällig weil ** (1/2) <=> ** 0 bzw. ** (3/2) <=> ** 1. Bei Dir funktionierte es aber nicht mehr, da Du korrekt r ** 3 berechnet hast, im Gravitationsgesetz muss es aber natürlich im Nenner r ** 2 sein.
Edit: Kommando zurück. In der vektoriellen Form muss natürlich doch r ** 3 im Nenner stehen. Die bisher schönen Grafiken dürften deshalb wohl einfach falsch schön sein. Man muss wohl einfach ein paar andere Planetendaten ausprobieren. Hier ein etwas besser funktionierendes Beispiel mit halbwegs attraktiven Grafiken: http://paste.pocoo.org/show/4022/. Die Skalierung der Kraftvektoren ist allerdings noch nicht ideal, so dass diese z.T. etwas wüst verlaufen.
MfG
HWK
Edit: Kommando zurück. In der vektoriellen Form muss natürlich doch r ** 3 im Nenner stehen. Die bisher schönen Grafiken dürften deshalb wohl einfach falsch schön sein. Man muss wohl einfach ein paar andere Planetendaten ausprobieren. Hier ein etwas besser funktionierendes Beispiel mit halbwegs attraktiven Grafiken: http://paste.pocoo.org/show/4022/. Die Skalierung der Kraftvektoren ist allerdings noch nicht ideal, so dass diese z.T. etwas wüst verlaufen.
MfG
HWK
So, jetzt ist auch das Problem mit der Vektordarstellung gelöst. Plot stellt scheinbar Vektoren in der Form (x y dx dy) dar und nicht (x1 y1 x2 y2). Ich habe eine Variante geschrieben, wo die Skalierung statisch im Programm festgelegt wird, http://paste.pocoo.org/show/4033/ und eine Version mit dynamischer Skalierung je nach der größten Kraft http://paste.pocoo.org/show/4032/. Aber auch diese Variante scheitert bei großen Kräften zwischen Punkten am Bildschirmrand.
MfG
HWK
Edit: Ich habe noch eine Variante mit Numpy geschrieben: http://paste.pocoo.org/show/4101/. Diese ist nicht so schön wie birkenfelds Vorschlag, aber auch nicht so schlecht. Einen Geschwindigkeitsgewinn kann ich subjektiv bisher nicht erkennen. Möglicherweise ist sie aber bei deutlich mehr Planeten effektiver.
MfG
HWK
Edit: Ich habe noch eine Variante mit Numpy geschrieben: http://paste.pocoo.org/show/4101/. Diese ist nicht so schön wie birkenfelds Vorschlag, aber auch nicht so schlecht. Einen Geschwindigkeitsgewinn kann ich subjektiv bisher nicht erkennen. Möglicherweise ist sie aber bei deutlich mehr Planeten effektiver.
Ein "physikalischer Fehler" entstand dadurch, dass die Vektoren berechnet wurden, nachdem a, v und s bereits geändert wurden. Diese waren deshalb nicht mit F konsistent. Es muss somit korrekt heißen:
Nach dieser Änderung sehen die Kraftvektoren wesentlich glaubhafter aus. Am besten sieht man das bei lediglich 2 Planeten.
MfG
HWK
Code: Alles auswählen
for i, p in enumerate(planeten):
# Die nächsten 4 Zeilen müssen vor den letzten 3 stehen
plot[i].append([p.s.x, p.s.y])
if iteration % 20 == 0:
vektor[i].append(
[p.s.x, p.s.y, F[i].x / 250, F[i].y / 250])
p.a = F[i] / p.m
p.v += p.a * dt
p.s += p.v * dt
MfG
HWK
-
- User
- Beiträge: 9
- Registriert: Dienstag 19. Dezember 2006, 21:13
Am rande bemerkt...
vielleicht kann decimal die fehlende genauigkeit der floats liefern. Bei Python2.5 unter Maxosx10.4 ist decimal wesentlich schneller als Numpy.
Andreas
vielleicht kann decimal die fehlende genauigkeit der floats liefern. Bei Python2.5 unter Maxosx10.4 ist decimal wesentlich schneller als Numpy.
Andreas