Newtons Gravitation mit GNUPlot und Python

Stellt hier eure Projekte vor.
Internetseiten, Skripte, und alles andere bzgl. Python.
Benutzeravatar
Agroschim
User
Beiträge: 28
Registriert: Sonntag 16. September 2007, 15:19

Newtons Gravitation mit GNUPlot und Python

Beitragvon Agroschim » Sonntag 16. September 2007, 17:19

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.
Benutzeravatar
Agroschim
User
Beiträge: 28
Registriert: Sonntag 16. September 2007, 15:19

Beitragvon Agroschim » Sonntag 16. September 2007, 17:21

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

[Code ausgelagert]

Edit (Leonidas): Code ausgelagert.
Benutzeravatar
birkenfeld
Python-Forum Veteran
Beiträge: 1603
Registriert: Montag 20. März 2006, 15:29
Wohnort: Die aufstrebende Universitätsstadt bei München

Beitragvon birkenfeld » Sonntag 16. September 2007, 18:55

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.
Dann lieber noch Vim 7 als Windows 7.

http://pythonic.pocoo.org/
Benutzeravatar
Agroschim
User
Beiträge: 28
Registriert: Sonntag 16. September 2007, 15:19

Beitragvon Agroschim » Sonntag 16. September 2007, 20:46

Ah... So geht das... :D
Benutzeravatar
HWK
User
Beiträge: 1295
Registriert: Mittwoch 7. Juni 2006, 20:44

Beitragvon HWK » Sonntag 16. September 2007, 21:50

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
Benutzeravatar
Agroschim
User
Beiträge: 28
Registriert: Sonntag 16. September 2007, 15:19

Beitragvon Agroschim » Montag 17. September 2007, 14:46

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.


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/
Benutzeravatar
birkenfeld
Python-Forum Veteran
Beiträge: 1603
Registriert: Montag 20. März 2006, 15:29
Wohnort: Die aufstrebende Universitätsstadt bei München

Beitragvon birkenfeld » Montag 17. September 2007, 14:55

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.


Das ist vollkommen in Ordnung. Beide Planeten wechselwirken und spüren beide eine Kraft (siehe Newtons zweites Gesetz).
Dann lieber noch Vim 7 als Windows 7.

http://pythonic.pocoo.org/
Benutzeravatar
birkenfeld
Python-Forum Veteran
Beiträge: 1603
Registriert: Montag 20. März 2006, 15:29
Wohnort: Die aufstrebende Universitätsstadt bei München

Beitragvon birkenfeld » Montag 17. September 2007, 14:56

Agroschim 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... :(


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
Dann lieber noch Vim 7 als Windows 7.

http://pythonic.pocoo.org/
Benutzeravatar
Agroschim
User
Beiträge: 28
Registriert: Sonntag 16. September 2007, 15:19

Beitragvon Agroschim » Montag 17. September 2007, 17:16

*bang* ! Ein Keks für den Herren dort! Danke, ich fall' vor die auf die Knie!
Benutzeravatar
Agroschim
User
Beiträge: 28
Registriert: Sonntag 16. September 2007, 15:19

Beitragvon Agroschim » Montag 17. September 2007, 18:01

Probleme gelöst, hier ist der Code:

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

Und hier ein Beispiel:

http://freenet-homepage.de/agroschim/planeten.png
Benutzeravatar
HWK
User
Beiträge: 1295
Registriert: Mittwoch 7. Juni 2006, 20:44

Beitragvon HWK » Montag 17. September 2007, 20:14

birkenfeld hat geschrieben: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
Zuletzt geändert von HWK am Montag 17. September 2007, 20:44, insgesamt 1-mal geändert.
Benutzeravatar
Agroschim
User
Beiträge: 28
Registriert: Sonntag 16. September 2007, 15:19

Beitragvon Agroschim » Montag 17. September 2007, 20:40

Interessant, könntest Du mir den Code ausdröseln? Wie gesagt, bin Python Neuling und generell nicht so der Programmierer...
Benutzeravatar
HWK
User
Beiträge: 1295
Registriert: Mittwoch 7. Juni 2006, 20:44

Beitragvon HWK » Montag 17. September 2007, 21:00

[code=] 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]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
Benutzeravatar
birkenfeld
Python-Forum Veteran
Beiträge: 1603
Registriert: Montag 20. März 2006, 15:29
Wohnort: Die aufstrebende Universitätsstadt bei München

Beitragvon birkenfeld » Montag 17. September 2007, 21:24

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)
Zuletzt geändert von birkenfeld am Dienstag 18. September 2007, 08:03, insgesamt 1-mal geändert.
Dann lieber noch Vim 7 als Windows 7.

http://pythonic.pocoo.org/
Benutzeravatar
HWK
User
Beiträge: 1295
Registriert: Mittwoch 7. Juni 2006, 20:44

Beitragvon HWK » Montag 17. September 2007, 21:31

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

Wer ist online?

Mitglieder in diesem Forum: 0 Mitglieder