Newtons Gravitation mit GNUPlot und Python

[Agroschim]

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.

[Agroschim]

Und hier, der Übersichthalber in einem zweiten Post, die zweite Variante:

[Code ausgelagert]

Edit (Leonidas): Code ausgelagert.

[birkenfeld]

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.

[Agroschim]

Ah... So geht das...

[HWK]

Die ganzen

Code: Alles auswählen

for i in range(anzplan):

sind sehr unpythonisch. Besser wäre z.B. für

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

oder

Code: Alles auswählen

for f_ in f:
    f_[:] = [0, 0]

Die while True-Schleife könnte man besser als

Code: Alles auswählen

for count in range(interationen):

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

[Agroschim]

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.

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...

http://paste.pocoo.org/show/3986/

[birkenfeld]

Die Kraft zwischen 2 Planeten A und B wird doppelt berücksichtigt, einmal von A nach B, dann von B nach A.

Das ist vollkommen in Ordnung. Beide Planeten wechselwirken und spüren beide eine Kraft (siehe Newtons zweites Gesetz).

[birkenfeld]

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...

Möglicherweise liegt das daran, dass

Code: Alles auswählen

r[i]**1/2

nicht das tut, was es soll, da "**" eine höhere Priorität als "/" hat.

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

[Agroschim]

*bang* ! Ein Keks für den Herren dort! Danke, ich fall' vor die auf die Knie!

[Agroschim]

Probleme gelöst, hier ist der Code:

http://paste.pocoo.org/show/3994/

Und hier ein Beispiel:

http://freenet-homepage.de/agroschim/planeten.png

[HWK]

Das ist vollkommen in Ordnung. Beide Planeten wechselwirken und spüren beide eine Kraft (siehe Newtons zweites Gesetz).

Stimmt. Da die Kräfte für beide Planeten aber bis auf das Vorzeichen identisch sind, kann man dann die Anzahl der Berechnungen halbieren:

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

Das sollte insbesondere bei vielen Planeten doch zu einer Zeitersparnis führen.
Für r und ds sind übrigens keine Listen notwendig.
MfG
HWK

[Agroschim]

Interessant, könntest Du mir den Code ausdröseln? Wie gesagt, bin Python Neuling und generell nicht so der Programmierer...

[HWK]

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

Du hast alle Kräfte außer den mit x bezeichneten bestimmt. Ich berechne nur die eine Hälfte (für i) und ordne die negative Kraft dem anderen Planeten (ii) zu. Der Rest entspricht Deinem Code, ist nur etwas optimiert: r, ds0, ds1 waren bei Dir unnötigerweise Listen. Das alte

Code: Alles auswählen

fg / r[i] ** (1 / 2)

kann man noch weiter zusammenfassen zu

Code: Alles auswählen

fg = g * m[i] * m[ii] / r ** (3 / 2)

Am Rande:

Code: Alles auswählen

f[i][0] += temp

entspricht wie in C

Code: Alles auswählen

f[i][0] = f[i][0] + temp

Code: Alles auswählen

ds0 * ds0

dürfte schneller sein als

Code: Alles auswählen

ds0 ** 2

Die kleinen Unterschiede in den Grafen dürften auf Rundungsfehlern bei der Integerdivision beruhen. Korrekter wäre es sicherlich, komplett mit Floats zu rechnen.
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()

wäre

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()

Hierbei taucht dann am Dateiende auch kein ']]' auf wie bei Dir. Entsprechendes gilt natürlich auch für das Schreiben der Vektor-Dateien.
MfG
HWK

[birkenfeld]

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)

[HWK]

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:

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 += ", "

MfG
HWK

[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

[Agroschim]

Hui, hier wird einem ja geholfen! Ich arbeite mich die Tage mal dadurch, heute finde ich keine Zeit mehr dafür! Vielen Dank jedenfalls schon mal!

[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.

[HWK]

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:

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

Nach dieser Änderung sehen die Kraftvektoren wesentlich glaubhafter aus. Am besten sieht man das bei lediglich 2 Planeten.
MfG
HWK

[andreas1965]

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