PDA

Archiv verlassen und diese Seite im Standarddesign anzeigen : Sonnensystem: Position der Planeten?



SRMeister
01.04.2011, 21:10
Huhu liebe User,
ich habe mir ein Programm geschrieben, welches die Planetenbahnen des Sonnensystems simuliert.
Nun bin ich soweit dass es recht gut funktioniert, aber leider fehlen mir informationen bezüglich der Positionen der Planeten. Bei Wiki bekommt man nur Abstand, Masse, usw. also alle wichtigen Parameter nur wo bekomm ich zusätzlich die aktuellen Positionen her? Um die Bahn korrekt festzulegen bräuchte man denke ich zusätzlich die Geschwindigkeit der Planeten?

Wenn bereits Interesse besteht könnte ich das Programm schon veröffentlichen, zb falls jemand helfen möchte, im Moment müssen die meisten Änderungen halt noch im Quellcode vorgenommen werden.(Man kann als Benutzer also noch nicht allzuviel machen) Ist geschrieben in C++.Net

Ziel soll u.a. sein, kurzfristig Asteroidenbahnen zu simulieren, bspw. Apophis, also bräuchte man auch deren Bahnparameter...

Stefan

Bernhard
01.04.2011, 22:24
Hi Stefan,

bei diesem Thema wirst Du Dich mit Programmen wie XEphem (Linux), Cartes Du Ciel und Stellarium (habe ich noch nicht ausprobiert) messen müssen. Sollte Dein Programm dabei durchfallen, empfehle ich das Buch von Oliver Montenbruck: "Grundlagen der Ephemeridenrechnung". Dort finden sich die wichtigsten Grundlagen und auch die Bahnparameter der Planeten.

Ob man Apophis ohne aufwändige Störungsrechnungen angehen sollte, wage ich zu bezweifeln.
Gruß

Kibo
01.04.2011, 22:36
Zum Thema aktuelle Position hab ich hier eine kleine Zusammenstellung der Onlineplanetarien (http://astro.nineplanets.org/astrosoftware.html#www)

Ich
01.04.2011, 22:49
Von Kibos Links empfehle ich Horizons (http://ssd.jpl.nasa.gov/?horizons). Das sind sozusagen die "offiziellen" Positionen.

SRMeister
02.04.2011, 00:07
Hi Bernhard, Kibo und Ich,

erstmal Danke für die hilfreichen Antworten.
@Bernhard:
Ich betreibe keine Ephemeridenrechnung. Ursprünglich habe ich ein Programm geschrieben zur N-Body Simulation von Protoplanetensystemen. Ich vermute einfach mal das ist eine andere Herangehensweise als Ephemeridenrechnungen.



Ob man Apophis ohne aufwändige Störungsrechnungen angehen sollte, wage ich zu bezweifeln.
Gruß
Störungsrechnungen sind natürlich da mit inbegriffen da jeder Körper mit jedem Anderen wechselwirkt und kollidieren kann.

Leider musste ich feststellen dass selbst mein Quadcore 3,2Ghz nicht annährend ausreicht um auch nur wenige Körper über lange Zeiträume zu simulieren. 10.000 Jahre für ca 100 Körper wäre beispielsweise ausreichend um diverse interessante Sachen zu simulieren, die so in jungen Planetensystemen abgehen. Aber leider komme ich selbst mit Multithreading also perfekter Ausnutzung aller 4 Kerne nicht annährend da ran :(

Allerdings für Zeiträume<1000Jahre und N<20 oder so, da lässt sich was starten.
Übrigens bei guter Genauigkeit ich denke Zeitschritte von 1-10sec sind machbar.

So jetzt schau ich erstmal wie ich die Planetenpositionen rein kriege.

mac
02.04.2011, 00:29
Hallo SRMeister,

UMa hatte mir vor einiger Zeit ein Verfahren skizziert, mit dem man solche Berechnungen erheblich beschleunigen und/oder genauer machen kann: http://www.astronews.com/forum/showthread.php?p=35493#post35493 ff

http://de.wikipedia.org/wiki/Compute_Unified_Device_Architecture könnte für Dich auch interessant sein.

Herzliche Grüße

MAC

SRMeister
02.04.2011, 01:20
Hallo MAC,
ich habe das extra in C++.Net geschrieben um später den "Kern" des Programms in C++ zu schreiben und damit den .Net Vorteil der einfachen Windows Programmierung mit dem Geschwindigkeitsvorteil von C++ zu kombinieren.
Insofern hab ich mir auch schon viele Gedanken über CUDA gemacht, was in .Net nicht machbar ist. Immerhin scheint das ja erheblich schneller als die CPU zu sein. Aber das bedeutet natürlich wieder extremen Aufwand die neue Art zu Programmieren zu lernen. Also mal schauen!
Die von UMa gemachten Vorschläge habe ich jetzt ehrlich gesagt auf Anhieb nicht verstanden. Allerdings habe ich wohl soviel verstanden, dass er letztendlich variable Zeitschritte, abhängig vom erwarteten Fehler vorschlägt?
Ich hatte mir diesbezüglich auch bereits Gedanken gemacht allerdings eher in der Form Zeitschritt abhängig von Entfernung zum nächsten Körper. So können nahe Körper mit hoher Genauigkeit berechnet werden... was im Endeffekt auf dasselbe hinausläuft !? Muss mir dazu noch mehr Gedanken machen :)
Jedenfalls bin ich sehr zufrieden (genau wie du damals :) ) dass es erstmal anständig läuft und für "einfache" Sachen wie Asteroidenbahnen sollte es erstmal so zu gebrauchen sein. Noch eine Frage, du schriebst, du benutzt den halben Zeitschritt ebenfalls, wie von UMa vorgeschlagen, aber in anderem Zusammenhang. Ich benutze nirgends halbe Zeitschritte. Wozu dienen die? War das nur auf die Ermittlung des Fehlers bezogen?

@Bernhard:
Habe mich eben noch etwas mit Bahnstörungen (http://de.wikipedia.org/wiki/Bahnst%C3%B6rung)beschäftigt wie von Dir vorgeschlagen, aber sehe erstmal keinen Anlass andere nicht-gravitative Effekte mit reinzunehmen? Es spielt sich alles auch weit genug von der Sonne ab so dass relativistische Effekte wahrscheinlich auch vernachlässigbar sind?

Ups eine Sache hab ich noch vergessen:
In dem von euch vorgeschlagenen HORIZONS Katalog(Danke nochmal, genau was ich gesucht habe) werden perfekterweise die Objekte als Position und Geschwindigkeit angegeben.
Hier als Beispiel die Erde, Koordinatenursprung ist die Sonne, Einheit ist AU:
JDCT
X Y Z <-- Position
VX VY VZ <-- Geschwindigkeit
LT RG RR LT = Lichtzeit(day); RG = Distance from Center(AU); RR = Range Rate(Änderungsrate von RG pro Tag (AU/day)

2455652.500000000 = A.D. 2011-Apr-01 00:00:00.0000 (CT)
-9.855554845305842E-01 -1.864518242419004E-01 1.975847807252305E-05 <-- Position
2.933461796767275E-03 -1.697932543422671E-02 2.779216797989756E-07 <-- Geschwindigkeit
5.793060517472141E-03 1.003037335417702E+00 2.739049035701513E-04

Einfach genial und genial einfach!


Grüße und gute N8 nun :)

Bernhard
02.04.2011, 07:59
@Bernhard:
Habe mich eben noch etwas mit Bahnstörungen (http://de.wikipedia.org/wiki/Bahnst%C3%B6rung)beschäftigt wie von Dir vorgeschlagen, aber sehe erstmal keinen Anlass andere nicht-gravitative Effekte mit reinzunehmen? Es spielt sich alles auch weit genug von der Sonne ab so dass relativistische Effekte wahrscheinlich auch vernachlässigbar sind?
Hi Stefan,

mit den Bahnstörungen meinte ich eher die gravitativen Einflüsse der Planeten untereinander. Man kann diese Einflüsse für vergleichsweise einfache Rechnungen auch vernachlässigen und bekommt damit Programme, die für rein visuelle Beobachter erst mal ausreichen. Es geht dabei aber mehr um die Frage, wo man einen Planeten abends am Sternhimmel sehen kann.
Gruß

SRMeister
02.04.2011, 13:51
So habe ich zuerst auch gedacht, habs mir aber trotzdem durchgelesen und naja, ich bin mir etwas unsicher ob relativistische Einflüsse vernachlässigt werden können, denn dort steht ja immerhin bei Merkur sind sie messbar.
Bei Asteroiden könnte auch der Sonnenwind eine Rolle spielen.

Ist also erstmal noch nix entschieden.

Bynaus
02.04.2011, 16:06
Bei Asteroiden könnte auch der Sonnenwind eine Rolle spielen.

Der Sonnenwind nicht, aber die Strahlung von der Sonne schon. Beim Yarkovsky-Effekt wandern die Asteroiden (bis zu einer Grösse von maximal etwa 30, 40 km, darüber ist das vernachlässigbar) aufgrund ihres Rotationsverhaltens von ihren anfänglichen Bahnen weg. Der YORP-Effekt hingegen beschleunigt die Rotation der Asteoriden und hat somit einen indirekten Einfluss auf den Yarkovsky-Effekt. Es gibt dazu ein nettes Review-Paper von Bottke et al, 2006, das bei "Annual Reviews in Earth & Planetary Science" erschienen ist.

Für ganz kleine Partikel (Mikrometeoriten) gibt es den Poynting-Robertson-Effekt, und für noch winzigere Partikel, die eine nennenswerte Ladung tragen können, kann auch der Sonnenwind eine Rolle spielen ("Blow-out").

mac
03.04.2011, 01:50
Hallo SRMeister,

ja, UMa steckt entweder so tief in der Materie, daß er vergessen hat wo für Anfänger die Stolpersteine rumliegen, oder er überschätzt meine Fähigkeiten, oder will mich dazu bringen meinen eigenen Grips zu bemühen. Letzteres ist ihm durchaus gelungen und im Nachhinein war das sicher nicht der schlechteste Weg. :)



Voraussetzung ist, daß Dein Programm richtig rechnet! Wie man das testet, beschreibe ich weiter unten - ist vielleicht nicht nötig, aber Deine Frage nach dem halben Zeitschritt, hat mich da in meiner Einschätzung dazu etwas verunsichert, weil Du eigentlich über dieses Problem gestolpert sein mußt.


Doch zunächst zu der von UMa beschriebenen Methode:
Du berechnest die Neue Position des Himmelskörpers A genau so, wie Du es sonst (hoffentlich richtig) machst. Sagen wir mal 1 Tag später. Du wiederholst exakt diese Berechnung, dieses mal aber mit kürzeren Zeitschritten, sagen wir mit 12 Stunden, dann brauchst Du schon zwei Zeitschritte bis Du bei der Position nach einem Tag angekommen bist und Dein Ergebnis für Position und Geschwindigkeit ist genauer. Du machst das nochmal, wieder für den selben Zeitabschnitt dieses mal aber unterteilt in 8 Stunden Schritte, also drei Schritte. Und so weiter (Schrittweite nach Fibonacci verändern ist da wahrscheinlich am besten, muß aber nicht, es geht auch mit anderen Verkürzungsregeln)

Wenn Du Dir nun die von Dir errechneten Werte für die X-Position anschaust und sie gegen die Schrittzahl aufzeichnest wirst Du sehen, daß diese X-Werte gegen einen bestimmten X-Wert konvergieren. Irgendwann werden Deine Ergebnisse, auf welcher X-Position der Planet nach einem Tag steht nicht mehr besser werden, sich nicht mehr verändern, egal wie klein Du Deine Zeitschritte machst.

Das Selbe gilt für die Y und Z Position und auch für die Geschwindigkeitsanteile in X, Y und Z Richtung.

Du kannst nun aus dem Kurvenverlauf, den diese X-Werte gegen die Schrittzahl zeigen ermitteln, wo dieser Konvergenzwert liegt/liegen wird, wenn Du die Schrittzahl weiter vergrößerst (also die Schrittweite weiter verkleinerst) ohne daß Du das bis zu eben dem oben beschriebenen Exzess wirklich durchrechnen mußt.

Und das genau ist der ‚Trick‘.

Ich selber habe das bisher nur mit wenigen Rechnungen ausprobiert und mich noch nicht darum kümmern können, ob ich wirklich in jedem Einzelfall eine eigene Kurvenanpassung rechnen muß, oder ob es ein ganz einfaches Verfahren gibt, den Konvergenzwert zu ermitteln. Ersteres kann ich zwar (war ein Werkzeug in meiner Diplomarbeit), aber es ist ziemlich aufwändig es auch zu programmieren (was ich müßte, wenn ich es in mein Programm einbauen will). Für Letzteres fehlen mir entweder die Erinnerungen wie man das macht, oder ich hab‘ es nie gelernt und selber bin ich bisher nicht (erneut?) drauf gekommen.



Nun zur halben Schrittweite.

Nehmen wir an, Du startest mit der Erde bei Z=0, Y=1AE, X=0, Deine Zeitschrittweite soll 1 Tag sein, dann gibt es mehrere Möglichkeiten wie Du die neue Position und Geschwindigkeit der Erde, nach diesem Tag ausrechnen kannst.

1. Du rechnest erst den neuen Geschwindigkeitsvektor (nach einem Tag Beschleunigung) mit dem Beschleunigungsvektor an dieser Position aus, und läßt die Erde mit diesem neuen Geschwindigkeitsvektor einen Tag lang Strecke machen. Das ist falsch.

2. Du läßt sie mit dem derzeitigen Geschwindigkeitsvektor einen Tag lang Strecke machen und rechnest an der nun nach einem Tag erreichten Position mit dem hier geltenden Beschleunigungsvektor den neuen Geschwindigkeitsvektor aus. Das ist auch falsch!

3. Du rechnest den neuen Geschwindigkeitsvektor mit dem hiesigen Beschleunigungsvektor aus, läßt sie aber einen halben Tag lang mit dem alten Geschwindigkeitsvektor bewegen und den zweiten halben Tag mit dem neuen Geschwindigkeitsvektor. (Man kann das nach diesem Startschritt bei den folgenden Schritten auch geeignet zusammenfassen)

Der für mich zuverlässigste Test ist: Sonne bleibt von jeder Beschleunigung verschont und die Erde bewegt sich auf einer Kreisbahn mit einer Geschwindigkeit von

V = Wurzel(G*mSonne/AE) = erste kosmische Geschwindigkeit

Nach jedem Schritt (hier ein Tag) muß ihr Abstand wieder 1AE sein.

Laß sie ruhig mal mit verschiedenen Schrittweiten einige tausend mal kreisen und schau Dir an, was passiert.

Herzliche Grüße

MAC

ispom
06.04.2011, 12:41
interessant finde ich auch, daß es Asteroiden gibt, die eine hufeisenförmige Bahn beschreiben...


Earth shares its orbit around the Sun with an asteroid in an exotic horseshoe-shaped orbit, say astronomers....
Apostolos Christou and David Asher at the Armagh Observatory in Northern Ireland say they've found one--an asteroid called 2010 SO16

http://www.technologyreview.com/blog/arxiv/26608/

Runzelrübe
06.04.2011, 13:23
Voraussetzung ist, daß Dein Programm richtig rechnet! Wie man das testet, beschreibe ich weiter unten

[...]

Der für mich zuverlässigste Test ist: Sonne bleibt von jeder Beschleunigung verschont und die Erde bewegt sich auf einer Kreisbahn [...]

Nach jedem Schritt (hier ein Tag) muß ihr Abstand wieder 1AE sein.


Mal aus reiner Neugier (denn ich entwickle Simulationsumgebungen beruflich):

Ihr nehmt jetzt aber hoffentlich nicht wirklich an, dass ihr in einer kontinuierlichen Simulation, in der ihr Ergebnisse als Eingangswerte weiterreicht, valide Werte mittels Double Precision rausbekommt, oder?! Bitte stellt sicher, dass ihr Datentypen verwendet, die sich um ihre Fehlerfortpflanzung kümmern. Informiert euch dazu bitte über robust geometric computing (http://www.cs.nyu.edu/yap/bks/egc/). Das allerschlimmste an Simulationen ist das Vertrauen in den Computer! :eek:

mac
06.04.2011, 16:22
Hallo Runzelrübe,


Ihr nehmt jetzt aber hoffentlich nicht wirklich an, dass ihr in einer kontinuierlichen Simulation, in der ihr Ergebnisse als Eingangswerte weiterreicht, valide Werte mittels Double Precision rausbekommt, oder?! Nein.



Bitte stellt sicher, dass ihr Datentypen verwendet, die sich um ihre Fehlerfortpflanzung kümmern.warum? Wird das Ergebnis dann genauer? Oder nehmen mir die nur die Fehlerrechnungen ab?



Informiert euch dazu bitte über robust geometric computing (http://www.cs.nyu.edu/yap/bks/egc/).danke für den Link!



Das allerschlimmste an Simulationen ist das Vertrauen in den Computer! :eek:Ja!

Herzliche Grüße

MAC

Runzelrübe
06.04.2011, 17:13
Hi MAC,


warum? Wird das Ergebnis dann genauer?

Das ist sehr stark abhängig vom Modell und den verschiedenen Skalen der Parameter. Allerdings werden Deine Werte durch robuste Datentypen immer auch für das mathematische Modell gültig. :cool:

Gerade bei numerischen Simulationen, die keinen sofort berechenbaren Zwischenzustand zu einem x-beliebigen Zeitpunkt erlauben, ist man so auf der sichereren Seite.

Du findest weitere Links, auch mit Negativbeispielen, auf: http://www.cs.nyu.edu/exact/resource/


The shear stresses were underestimated by 47%, leading to insufficient design. In particular, certain concrete walls were not thick enough



The effect of this inaccuracy on the range gate's calculation is directly proportional to the target's velocity and the length of the the system has been running. Consequently, performing the conversion after the Patriot has been running continuously for extended periods causes the range gate to shift away from the center of the target, making it less likely that the target, in this case a Scud, will be successfully intercepted.

In Simulationen, in denen es um Systeme geht, die ihre Einflüsse zu einem Vektor addieren, können wir nach zig Zehntausend Iterationen eine nette Fehlersumme zusammen haben. Jeder Körper einzeln führt seinen vorherigen Fehleranteil weiter in die jeweils falsche Richtung. Sieht stabil aus auf einer gewissen Skala und für einen gewissen Zeitraum, ist aber ungenügend für Aussagen, das Ende betreffend. Gerade das Beispiel mit der Scud-Rakete, durch die Menschen gestorben sind, sollte zum Nachdenken anregen.

MGZ
06.04.2011, 17:23
Das allerschlimmste an Simulationen ist das Vertrauen in den Computer! :eek:

Es gibt Methoden, um Simulationen sehr gruendlich zu testen. Vertrauen allein reicht niemals beim Programmieren.

SRMeister
06.04.2011, 23:24
Nur um nochmal eine kurze Rückmeldung zu geben:

@mac das mit den halben Schritten war mir in der Tat nicht bewusst. Aber du hast recht, meine Sonne ist im Moment verankert. Es muss also was falsch sein. Die Planetenbahnen sehen aber gut aus.

Ich habe die Geschwindigkeit jetzt nochmal verdoppeln können: Da die berechnete Kraft F(a-b) zwischen (a) und (b) gleich der Kraft F(b-a) ist spart man sich mit wenig Aufwand genau die Hälfte aller Berechnungen.

Werde aber wohl doch in Richung Cuda gehen müssen und das erklärte Ziel ist N*log(N) statt N²
wie?
hier mehr (http://beltoforion.de/barnes_hut/barnes_hut_de.html#idIntro)

die besten Grüße in die Runde
Stefan

mac
07.04.2011, 01:02
Hallo Stefan,

hör auf Runzelrübe.

Eine solche Simulation des Sonnensystems ist eine nette Spielerei, bei der man mit (vor)Freude auf’s Ergebnis viel Neues, ganz unauffällig nebenbei lernen kann. Ich hab‘ sowas 1985 mit den Daten des Sonnensystems aus meiner alten Logarithmentafel noch aus der Schulzeit, zum ersten Mal gemacht und hatte keine Ahnung, wie ich an die echten Positionen ran kommen könnte. Welch ein Schlaraffenland steht uns da mit dem Internet zur freien Verfügung! Ich hab‘ ‚Nemesis‘ durchs Sonnensystem rauschen lassen und mich über die durcheinander gewürfelten Planetenbahnen gefreut. :D Aber besonders darüber, wie genau ein solcher ‚Stern‘ treffen muß, um wirklich Schaden anzurichten. :)

Und dann, bei einem erneuten, noch genaueren Treffer, raste auf einmal unser Mond mit ungeheurer Geschwindigkeit davon. Was war passiert? Ich hatte mit festen Zeitschritten gerechnet. Das geht ganz unauffällig gut, solange es keine großartigen Entfernungsunterschiede zwischen den einzelnen Zeitschritten gibt, aber wenn man einen Tag lang mit einer Beschleunigung rechnet, die so höchstens ein paar Minuten dauert, dann bekommt man sehr drastisch die Grenzen solcher Simulationen vor Augen geführt und damals für mich die Grenzen, wie weit ich solche ‚Spielchen‘ in der Mittagspause ausbauen konnte und wie weit nicht.

Du kannst Deinen Simulator nach den Grundtests mit den echten Daten des Sonnensystems testen (aus dem Link von Ich. Laß Dir die Positionen und Geschwindigkeiten(auch die der Sonne)auf den Schwerpunkt des Sonnensystems bezogen ausgeben) und Du wirst auch sehen, daß Du allein mit Newton nicht auskommst, wenn es sehr genau sein soll und das über mehrere Jahre, besonders bei den inneren Planeten.

Schau Dir an, mit welcher Genauigkeit Du an die Massen der Himmelskörper kommst. Das sind alles gerundete Werte. Zum Glück sind die Umlaufbahnen mehr oder minder Kreisförmig, so daß sich die Fehler duch die Rechengenauigkeit bei den Vektor-Berechnungen nicht vollständig addieren. Wenn Du mit floating Point Daten arbeitest, dann sortiere Deine Ergebnisse vor den Additionen nach (absoluter, also vorzeichenloser) Größe, und fang die Additionen mit den kleinen Zahlen an, sonst fällt unnötig viel unter den Tisch. Hier wirst Du dann auch sehen, wo die Beschränkungen in der Rechengenauigkeit liegen und Du kannst Dir in erster Näherung einfach überschlagen, nach wieviel Iterationen Dein Rechenfehler wie groß geworden sein muß, und eben auch, daß zigtausende von Jahren nur noch Hausnummern liefern können.

Dein Plan es mit Barnes und Hut zu versuchen, ist gut, aber er vergrößert auch den systematischen Fehler - ok der tritt allerdings bei einer solchen Dynamik wie sie im Sonnensystem herrscht, eher chaotisch verteilt auf, so daß Du hoffen darfst daß es ein noch brauchbarer Kompromiß ist. Vergleiche einfach die Ergebnisse mit denen ohne diese Zusammenfassungen.

Herzliche Grüße

MAC

SRMeister
09.04.2011, 16:50
und Du wirst auch sehen, daß Du allein mit Newton nicht auskommst, wenn es sehr genau sein soll und das über mehrere Jahre, besonders bei den inneren Planeten.
Die Frage hab ich ja schon weiter vorne aufgewurfen.
Mein Projekt habe ich in Gedanken in zwei Teilprojekte unterteilt. Das Erste war eine möglichst langfrisitge bezogen auf Planetesimale und die Entstehung von Planeten daraus. Das ist vorerst gescheitert da ich noch kein Barnes Hut Algo (also logarithmische Skalierung) implementiert habe und auch noch nicht CUDA verwende.
Das Zweite Teilprojekt ist die Simulation des aktuellen Sonnensystems, um möglichst genaue Voraussagen treffen zu können, was mit einigen Asteroiden passiert. Dazu reichen Simulationszeiten von max. ca. 100 Jahren, was aber möglichst genau sein soll.
zur Genauigkeit: Ich rechne ausschließlich mit double, also 64bit. Double hat ca. 17 Nachkommastellen.

Schau Dir an, mit welcher Genauigkeit Du an die Massen der Himmelskörper kommst. Das sind alles gerundete Werte.
Ich denke in der Horizons Datenbank ist alles in ausreichender Genauigkeit zu finden.

Ich habe nun das Programm in "unmanaged C++" übersetzt weil sich .NET nicht mit CUDA C erweitern lässt und weil es wesentlich schneller läuft.
Bald gibts eine Veröffentlichung.

beste Grüße,
Stefan

Runzelrübe
09.04.2011, 19:49
Ich habe nun das Programm in "unmanaged C++" übersetzt weil sich .NET nicht mit CUDA C erweitern lässt und weil es wesentlich schneller läuft.

Ein Wort der Warnung: Sei bitte vorsichtig, dass Du damit keinen Schritt zurück machst, indem Du Dir durch eine schnellere Berechnung vorgaukelst, dass Dein Programm besser wird. Es wird eben nur eins: schneller. Jede Software kann man in Hardware gießen und genau das wird auf einer GPU für die auf dem Bildschirm benötigten Genauigkeitsvorgaben getan. GPUs arbeiten jedoch mit noch geringeren Bit-Breiten für Datentypen. Um es für Dich bildlich auszudrücken: ausgerechnet Deine kleinen Asteroiden werden immer ungenauer, weil ihre Wichtungsanteile unter den Tisch fallen.

Da ich ebenfalls in .NET entwickle, sei Dir noch ans Herz gelegt, Dich zum Thema unmanaged lieber in Richtung Marshalling zu informieren. Jede C/C++ Bibliothek kann durch DLLImport in managed Code integriert werden. Unmanaged ist schnell auch unsafe, daher heißt auch die Compilerdirektive genau so: unsafe. Sei bitte vorsichtig, sowohl mit Deinen Datentypen und auch mit Deiner Speicherverwaltung.

Worin ich Dich gern bestärken möchte: Plane Deine Schritte. Mach bitte einen Schritt nach dem Anderen. Mach kleine Schritte, diese dann allerdings vollständig. Und geh danach erst weiter bis zum Ende.

Ich hatte vor Jahren ebenfalls mal ein Programm (in C#.NET) geschrieben, welches mir .avi Videos zur Simulation von Sternbewegungen in Sternhaufen ausgespuckt hat und habe mir auch nicht reinreden lassen. Mittlerweile weiß ich es besser. Meine Ergebnisse waren nett anzuschauen, aber leider komplett unbrauchbar. Und das hat mich eher enttäuscht als motiviert. Ich kann Dir nur raten, red Dir nicht ein, Ratschläge selektiv zu ignorieren, weil Du das ja nur aus Hobby machst. Du könntest sonst am Ende ebenfalls enttäuscht sein, wenn Du das Programm präsentierst und es nur für 10-Körper-Probleme ähnlicher Größenordnungen gut läuft.

Bernhard
09.04.2011, 20:08
Worin ich Dich gern bestärken möchte: Plane Deine Schritte. Mach bitte einen Schritt nach dem Anderen. Mach kleine Schritte, diese dann allerdings vollständig. Und geh danach erst weiter bis zum Ende.
Hi Rübe,

man kann prinzipiell aber auch aus Fehlern lernen. Meine Devise lautet in diesem Fall deswegen: Lieber ein paar Fehler machen, als gar nix machen.

Dass .NET soviel besser ist, als Altbekanntes muss mir .NET erst mal beweisen.
Gruß

Runzelrübe
09.04.2011, 23:11
Hi Bernhard,


Dass .NET soviel besser ist, als Altbekanntes muss mir .NET erst mal beweisen.

Erklär mir bitte: hattest Du wirklich vor, das Thema mit diesem Kommentar jetzt in eine bestimmte Richtung (die arg nach einer Diskussion PC vs. MAC aussieht) zu lenken? Was meinst Du denn mit altbekannt? Und wieso ist das Altbekannte denn so viel besser? Weil Du es kennst? Ja, der Mensch ist schon ein faules Tier. Was macht das Altbekannte denn besser? Kleiner Tipp: altbekannt bedeutet zuerst mal, dass es früher dagewesen ist, nicht automatisch, dass es das beste ist, was es je geben wird. Es gab auch eine Evolution der Computersprachen. Die Erde war lange Zeit eine Scheibe, das war auch altbekannt, hatte aber ein paar Makel. Dann war sie kugelförmig, dann doch eher ein Geoid. Hey, hat sich doch gelohnt. Niemand wird mehr daherkommen und behaupten, sie wäre eckig. Und uns allerliebster Newton war auch zuerst da, doch dann kam Einstein. :rolleyes:

Daher: denkt bitte nicht, dass irgendeine x-beliebige Programmiersprache einer anderen überlegen ist, denn sie ist auch nur eine weitere x-beliebige Krücke zur Binärwelt. Ja komm, ich kenn euch: es ist so schön bequem, sich zurückzulehnen, wenn man eine Sprache endlich beherrscht. :D Manche dieser Binärkrücken sind sicher, manche hübsch, manche einfach, manche für Anfänger, manche nicht kompatibel zum Problem usw.

Aber: wenn ich lese, dass ein Programm, das in Hochsprache A geschrieben wurde, in Hochsprache B konvertiert wird, damit dies und das funktioniert, obwohl dies (fremd-)wissentlich nicht notwendig wäre, dann klingeln bei mir die Alarmglocken und ich schreite zur Tat.

Ich stimme Dir im Übrigen im Grundsatz zu, den Ansporn zu haben, aus Fehlern zu lernen. Das ist meiner Meinung nach fundamental richtig. Doch dazu gehört eben auch, die Fehler zu erkennen oder darauf aufmerksam gemacht zu werden. Und da liegt der Hase im Pfeffer. Man hat es so schön einfach, wegzuschauen.

Fehler zu wiederholen, führt zu noch mehr Fehlern. Und ohne qualifizierte Tests wird das nix. Diese fangen aber leider ganz unten an, beim Urschleim, dem ollen blöden Rechner mit seinen ollen blöden Registern, die nur gesagt bekommen, was sie wann zu tun haben. Schieb mal das hierher und dann dreh Dich im Kreis. Toll!

Gegenfrage: Führst Du mathematische Korrekturterme ein, damit die Messwerte stimmen, damit sie Dir der erstbeste Kollege später um die Ohren haut? Und zwar obwohl Dir gesagt wurde, dass diese Korrekturterme ein no-go sind? Manchmal muss man eben zwischendurch nochmal sein Wissen aufbessern oder nachhaken. War das denn jetzt wirklich so schlimm, das zu erwähnen? Oder worin liegt da das Problem, Deiner Meinung nach? Weil ich nicht ja und Amen dazu sage, dass SRMeister sich Mühe gibt? Hey, ich teile mein Wissen freiwillig, weil ich ihm damit helfen möchte. Das ist verdammt nochmal viel wert!

LG,
Runzel

SRMeister
09.04.2011, 23:37
Hallo,


GPUs arbeiten jedoch mit noch geringeren Bit-Breiten für Datentypen.
soweit ich weis stimmt das zumindest im Fall Cuda nicht, wenn man die ersten Cudafähigen Baureihen ausklammert.

Zur Verbesserung der GPU-Computing-Fähigkeiten verfügen die Grafikprozessoren der „Fermi-Architektur“ als erste überhaupt über eine komplette Unterstützung von C++ und sind, genau wie die Radeon-HD-5000-Serie von AMD, mit dem IEEE-754-2008-Standard vollständig kompatibel.


Da ich ebenfalls in .NET entwickle, sei Dir noch ans Herz gelegt, Dich zum Thema unmanaged lieber in Richtung Marshalling zu informieren.
Genau so ist es natürlich geplant. Ich hab im Moment ne c++ Konsolenanwendung und ne .NET Win Anwendung und die sollen zusammengeführt werden. Leider mach ich das auch das erste mal deswegen haperts etwas mit der Umsetzung. Stichwort: Wrapper Klasse erstellen. Habe ja nicht gewusst dass es noch andere .net nerds hier gibt :)
PS wie schon gesagt würde mich über jede Hilfe freuen :)


Sei bitte vorsichtig, sowohl mit Deinen Datentypen und auch mit Deiner Speicherverwaltung.
immer!


Meine Ergebnisse waren nett anzuschauen, aber leider komplett unbrauchbar. Und das hat mich eher enttäuscht als motiviert.
Warum waren sie unbrauchbar?


Ich kann Dir nur raten, red Dir nicht ein, Ratschläge selektiv zu ignorieren, weil Du das ja nur aus Hobby machst. Du könntest sonst am Ende ebenfalls enttäuscht sein, wenn Du das Programm präsentierst und es nur für 10-Körper-Probleme ähnlicher Größenordnungen gut läuft.
Natürlich denke ich über jeden gemachten Ratschlag nach. Es gibt in meinem Kopf eine ToDo Liste wozu bspw. auch die Durcharbeitung des verlinkten Textes über Fehlerfortpflanzung gehört und die Beschäftigung mit dem Fehler sowieso. Nur muss ich erst sehen ob ich das Programm ansich überhaupt soweit bekomme, dass es überhaupt Ergebnisse liefern kann. Es kann ja auch sein, dass bei 10 Körpern das Programm so langsam wird dass ich aufgebe. Dann wäre alle jetzige Mühe zum Thema Fehler verschwendet.



Dass .NET soviel besser ist, als Altbekanntes muss mir .NET erst mal beweisen.
Ich hab auch mit MFC angefangen. Aber ehrlich gesagt würd ich jetzt nie wieder freiwillig zu MFC zurückkehren. Und C++ wird nurnoch gemacht wenns nötig ist, wie hier.
NET ist vielleicht nich allgemein besser, dafür aber auf jeden Fall, nimmt es einem sehr viel arbeit ab, gerade im Bereich Windows Programmierung, aber eigentlich auf allen Gebieten.

Stefan

FrankSpecht
09.04.2011, 23:44
Moin, Runzel,

hattest Du wirklich vor, das Thema mit diesem Kommentar jetzt in eine bestimmte Richtung (die arg nach einer Diskussion PC vs. MAC aussieht) zu lenken?
Ich glaube nicht. Ich denke, Bernhard geht es hauptsächlich um das .NET (http://de.wikipedia.org/wiki/.NET).
Ist halt Microsoft-spezifisch und nicht ohne weiteres kompatibel zu Standard-C (oder C++), so dass ein solches Programm ohne Aufwand nicht auch auf anderen Betriebssystemen lauffähig ist.

Außerdem geht es um die Ausführungsgeschwindigkeit, die ja, je nach Aufgabe, auch eine Rolle spielen könnte:

Damit entfernt sich die Softwareentwicklung, insbesondere die Anwendungsentwicklung, insgesamt von dem Augenmerk auf Ausführungsgeschwindigkeit (http://de.wikipedia.org/wiki/Leistung_%28Informatik%29), die angesichts einer immer schneller werdenden Rechengeschwindigkeit (http://de.wikipedia.org/wiki/Floating_Point_Operations_Per_Second) zunehmend in den Hintergrund rückt.

Bernhard
10.04.2011, 00:00
hattest Du wirklich vor, das Thema mit diesem Kommentar jetzt in eine bestimmte Richtung (die arg nach einer Diskussion PC vs. MAC aussieht) zu lenken?
Hi RR,

mir ging es eigentlich nur darum SRMeisters lobenswerte Motivation vor überzogenen Software- und Codierungsrichtlinien zu schützen.
Gruß

Bernhard
10.04.2011, 09:09
Ich hab auch mit MFC angefangen. Aber ehrlich gesagt würd ich jetzt nie wieder freiwillig zu MFC zurückkehren. Und C++ wird nurnoch gemacht wenns nötig ist, wie hier.
Hi Stefan,

vielen Dank für den Tip, aber um mal wieder auf die Physik und Numerik zurück zukommen: Mit welchem Verfahren machst, bzw. willst Du eigentlich die Integration implementieren und in welchen Koordinaten rechnest Du? Ich könnte mir vorstellen, dass man in Kugelkoordinaten schon mal weiter kommt, als in euklidischen Koordinaten.
Gruß

Bernhard
10.04.2011, 09:33
3. Du rechnest den neuen Geschwindigkeitsvektor mit dem hiesigen Beschleunigungsvektor aus, läßt sie aber einen halben Tag lang mit dem alten Geschwindigkeitsvektor bewegen und den zweiten halben Tag mit dem neuen Geschwindigkeitsvektor. (Man kann das nach diesem Startschritt bei den folgenden Schritten auch geeignet zusammenfassen)
Hi MAC/Stefan,

eine Beschreibung dieses Verfahrens findet man u.a. auch hier: http://en.wikipedia.org/wiki/Midpoint_method.
Gruß

Runzelrübe
10.04.2011, 11:27
Hi Stefan,

die Einschränkung auf CUDA ist ähnlich bindend wie eine Sprache, mit dem Unterschied, dass zwangsläufig auch noch NVidia Hardware vorausgesetzt wird. Und dass CUDA jetzt endlich auch Double Precision kann, das war mir neu. Aber auch NVidia schläft eben nicht. :)

Zu Deinen Bedenken bezüglich Geschwindigkeit. Auf meinem damaligen AMD Athlon 2100+ liefen sowohl die Berechnung als auch 3D-Darstellung eines 60-Körper-Systems wie ein Echtzeitvideo. Einzig für die Simulation von 2500-Körper-Systemen habe ich eine Videoauslagerung integriert und meinen Rechner nachts arbeiten lassen.


Warum waren sie unbrauchbar?

Weil mir damals einige Parameter als vernachlässigbar klein erschienen. Sonnenwind und Gezeitenkräfte hatten meines Erachtens viel zu geringe Auswirkungen. Ich hatte die Auswirkungen naher Vorübergänge einfach aus Bequemlichkeit nichtmal durchgerechnet. Dazu kommt, dass die Ausbreitungsgeschwindigkeit der Gravitation bei mir instantan erfolgte, statt mit Lichtgeschwindigkeit. Ich benötigte immer ein schweres schwarzes Loch im Zentrum zur Stabilisierung eines Kugelsternhaufens. Das seltene Ereignis einer Nahbegegnung eines Sterns mit diesem hatte dann regelmäßig das Ende der schön anzusehenden Simulation zur Folge. Denn auch die Annihilation eines Sterns (und sei es nur durch einen unelastischen Stoß) hatte ich weggelassen. Ganz am Ende habe ich mir einfach gesagt, dass es ja ganz nett aussieht und habe mich anderen Dingen gewidmet. Einige Jahre später wusste ich dann besser Bescheid, habe mich nur leider nie wieder rangesetzt. Bei einem Thema konnte ich mir jedoch 100%ig sicher sein. Dass die verwendeten Datentypen mich nicht im Stich lassen konnten. Such mal nach Jonathan Richard Shewchuk. Prioritäten liegen eben bei jeder Person woanders. Das macht uns so schön unterschiedlich. :)

Schönes Beispiel zum Nachdenken:

(333.75−a^2)*b^6+a^2*(11*a^2*b^2−121*b^4−2)+5.5*b^ 8+a/(2*b)

a = 77617
b = 33096

Ergebnisse:

32-bit: 1.172604
64-bit: 1.1726039400531786
128-bit: 1.1726039400531786318588349045201838

Wie auch immer, das richtige Ergebnis lautet: −0.8273960.

-----------

Hi Bernhard,


mir ging es eigentlich nur darum SRMeisters lobenswerte Motivation vor überzogenen Software- und Codierungsrichtlinien zu schützen.

Fremdschutz ist sinnvoll, wenn derjenige, den es zu schützen gilt, bevormundet werden möchte oder diese Aufgabe jemandem übertragen hat. Mir kam es aber zumindest so vor, als wisse Stefan bereits, dass ich nicht gegen ihn argumentiere, sondern auf Fehlerquellen hinweise. Sein Programm kann definitiv auch mit Double Precision Berechnungen gute Ergebnisse liefern. Und wie er bereits schrieb, möchte er sich die Links später noch zu Gemüte führen.

SRMeister
21.04.2011, 22:23
so es geht weiter.

Ich habe nun eine normale C++ Konsolenversion gemacht und ihr dürft sie gerne zerpflücken - aber hoffentlich auch beim Verbessern helfen :D

Im Moment läuft das ganze so:
Euler Methode (Keine Zwischenpunktberechnung), Zeitschritt eine Sekunde, innere Planeten bis Erde plus Mond,
im angehängten Beispiel werden etwa 2 Monate integriert. Die Position stimmt mit der von NASA folgendermaßen nicht überein :

X Abweichung 24.500km, Z Abweichung 1800km, Y irgendwo darunter...
Nur als Referenz die Erde hat sich in der Zwischenzeit "Luftlinie" folgende Strecke bewegt:
X 870 Mio km, Y 110 Mio km, Z 4 Mio km

Der Fehler ist also erstmal "recht gering", vielleicht durch die fehlenden äußeren Planeten oder durch die Ungenauigkeit der Methode.

Der nächste Schritt ist der Wechsel weg von der Euler Methode hin zu was genauerem.
Wahrscheinlich erstmal die Midpoint Methode.
Danach wird das Ganze als DLL kompiliert und wieder in .NET eingebunden damit die Grafikmöglichkeiten wieder da sind.

Was mich interessieren würde: Wenn ich die Genauigkeit steigere, bspw. durch wechsel Euler -> Midpoint, verdoppelt sich die Rechenzeit. Das lohnt sich ja nur wenn die Genauigkeit mehr steigt als Faktor 2 damit ich die Schrittweite größer machen kann. Es gibt ja dann noch mehr Methoden (Runge-Kutta-Methoden). Allerdings muss es ja irgendwo ne Grenze geben, wo die Genauigkeit sich zum Rechenaufwand sich nicht mehr verbessert ??

hier für jeden zum Nachschauen. (http://www.coding-facility.com/sr/PSystemC.rar) Inklusive EXE Datei für alle die zu faul sind sich Visual Studio Express (kostenlos) zu saugen, oder keine Ahnung vom Programmieren haben :)

schönen Abend!
Stefan

Nathan5111
22.04.2011, 02:45
Schönes Beispiel zum Nachdenken:

(333.75−a^2)*b^6+a^2*(11*a^2*b^2−121*b^4−2)+5.5*b^ 8+a/(2*b)

a = 77617
b = 33096

Vielen Dank für dieses wunderschöne Beispiel, ich habe es (erfolgreich!) bei einem 'hoffnungslosen Fall' eingesetzt, es ließ sich mit 'Bordmitteln' lösen.

Bernhard
22.04.2011, 10:52
Ich habe nun eine normale C++ Konsolenversion gemacht und ihr dürft sie gerne zerpflücken - aber hoffentlich auch beim Verbessern helfen :D
Hi Stefan,

ich hatte oben leider etwas Panik bekommen, dass du deine Ankündigung einer Veröffentlichung von Code wieder fallen läßt. Deswegen erst mal vielen Dank für den Code. Vorallem auch deswegen, weil ich viele Jahre mit Visual C++ 6.0 vor allem beruflich gearbeitet habe und neulich mit ziemlichem Schrecken festgestellt habe, dass Visual C++ 2010 doch sehr viele Änderungen mit sich bringt. Der Hinweis auf unmanaged code war deswegen für mich eine mittlere Erlösung von unüberwindlich erscheindenden Migrationsproblemen.

Zusammen mit Internet und diversen Büchern werde ich mich jetzt erst mal in die Problematik einarbeiten, alte Programme fit für Windows7 zu machen.

Zu Deinem Programm kann ich nur sagen, dass viele Wege nach Rom führen und ich kann nur sehr empfehlen bei vorhandener Zeit auch mehrere Wege auszuprobieren. Man kann dann Ergebnisse bezüglich der Genauigkeit vergleichen. Der Genauigkeit wird in der Praxis zumindest in meinem Umfeld nicht selten der Vorzug gegeben im Vergleich zur Rechenzeit.
Gruß

SRMeister
22.04.2011, 12:11
Hey Bernhard,
ja mit VS6 hab ich auch angefangen. Um ehrlich zu sein hab ich auch ewig damit gearbeitet, noch als es schon längst VS2003.NET gab.
Das Projekt jetzt lässt bei mir auch teilweise gute Erinnerungen an die alten Zeiten aufleben, hab lange nix in unmanaged C++ gemacht.
Das jetzige Projekt wurde mit VS2008 erstellt, allerdings kannst du ganz einfach mit jedem beliebigen C++ Compiler ein Projekt erstellen mithilfe der paar Quellcode Dateien. Einzig die 2 standart-includes in main.cpp müsste man evtl. anpassen.

hier noch eine kleine Geschwindigkeitsverbesserung:

in der Datei PSystem.cpp, die Funktion DoStep, dort der Inhalt der if-Anweisung kann durch folgendes ersetzt werden. Bringt etwa 30% speed:

PVektor *tmp = PVektor::Sub(target->pos, source->pos); // Beschleunigung berechnen
double dist = tmp->Len();
tmp->Mul(target->mass/(dist*dist*dist));
source->Accel->Add(tmp); // Beschleunigung hinzufügen
delete tmp;

(Endlich bringt die Code Funktion hier mal was :D)
Übrigens dauern bei mir die 2 Monate etwa 14 sec auf 2,4 Ghz Intel Core 2 Prozessor, ich finde das ist recht schnell und wenn ich bei der ungefähren Geschwindigkeit bleibe und die Genauigkeit weiter steigern kann, kann ich denke mit dem NASA Horizons hoffentlich mal mithalten :)
Ziel wie gesagt, höchste Genauigkeit, Integration über 100 Jahre, Berechnung des inneren Sonnensystems inkl. einiger NEAs.

CUDA verschiebe ich erstmal, ist mir noch etwas zu weit entfernt (im Sinne von ich hab einfach zuwenig Ahnung davon) aber SSE2 ist noch ne alternative. Mal schauen.
Erstmal möchte ich aber die Midpoint Methode haben! Hoffe da noch auf Unterstützung, Codebeispiele oder ähnliches!

Stefan

Bernhard
22.04.2011, 17:49
Erstmal möchte ich aber die Midpoint Methode haben! Hoffe da noch auf Unterstützung, Codebeispiele oder ähnliches!
hier mal ein kleiner Code-Review:

a) Tippfehler in main.cpp, Zeile 57: Starzeit -> Startzeit

b) PSystem.cpp, Zeile 86: Die Deklarationen
PVektor *tmp und
double dist würde ich aus sämtlichen Schleifen herausziehen und an den Anfang der Funktion stellen. Zeile 92 kann man dann weglassen, weil der Speicher für die Variable beim Verlassen der Funktion automatisch freigegeben wird. Mag sein, dass der Compiler das eh so übersetzt, aber damit wird der Code auch übersichtlicher und lesbarer.

Insgesamt würde ich auch aus Gründen der Rechengeschwindigkeit unter dem if also eher den folgenden Code verwenden:


tmp = PVektor::Sub(target->pos, source->pos); // Beschleunigung berechnen
dist = tmp->Len();
tmp->Mul(Gconst*target->mass/(dist*dist*dist));
source->Accel->Add(tmp); // Beschleunigung hinzufügen


Damit hat man dann die Beschleunigung in SI-Einheiten. Die Gravitationskonstante Gconst sollte IMO nur an dieser Stelle in die Rechnung einfließen. Zeile 94 würde ich rausnehmen und in die Funktion Planet::Move einbauen, damit die Variablennamen auch wirklich das bezeichnen, was sie speichern. D.h. die Variable accel sollte nur Beschleunigungen und keine Positionen speichern.

Der Funktion Move würde ich deshalb die Schrittweit als Parameter übergeben und den Code in der Funktion wie folgt anpassen:


Velocity->Add(Accel * schrittweite);
pos->Add(Velocity * schrittweite);
Accel->Reset();


Man kann dann das Projekt als Ganzes kopieren und in der Methode Move die Midpoint-Methode implementieren. Eventuell sollte man die Methode Move auch umbenennen. Wie wäre es z.B. mit Planet::GetNextPosition(double schrittweite)?

c) Bei der Klasse Planet fehlt mir noch der Destruktor, der den mit new reservierten Speicher für die Vektoren pos, Accel und Velocity wieder freigibt.
Gruß

SRMeister
22.04.2011, 21:01
Zeile 92 kann man dann weglassen, weil der Speicher für die Variable beim Verlassen der Funktion automatisch freigegeben wird. Mag sein, dass der Compiler das eh so übersetzt, aber damit wird der Code auch übersichtlicher und lesbarer.
Das ist in dem Fall nicht korrekt da es sich um Zeiger handelt und auch wenn die Zeigervariable selbst gelöscht wird so bleibt der Inhalt des Speichers erhalten, weshalb ein delete erforderlich ist.

Alle anderen vorgeschlagenen Sachen sind mMn korrekt und deshalb umgesetzt.


Die Gravitationskonstante Gconst sollte IMO nur an dieser Stelle in die Rechnung einfließen.
Aus Gründen der "mathematischen Correctness" wäre das richtig, aber da bei dieser Variante die Gconst mit multipliziert werden muss, hab ich sie gleich beim Festlegen der Masse einbezogen. So spart man sich die eine Multiplikation dort wo es auf speed ankommt.


Zeile 94 würde ich rausnehmen und in die Funktion Planet::Move einbauen, damit die Variablennamen auch wirklich das bezeichnen, was sie speichern. D.h. die Variable accel sollte nur Beschleunigungen und keine Positionen speichern.

source->Accel->Mul(TconstSquare); // Pos = Zeit*Zeit*Kraft

ist meiner Meinung nach korret. Folgendermaßen:
Ich habe ja weiter vorne erwähnt dass ich Quasi rausgefunden habe: Normalerweise berechnet man die Kraft auf einen Körper. Diese Kraft ist aber auf den anderen genau gleich groß, deshalb spart man sich an der Stelle eine Berechnung. F(hin) = F(zurück).
Das habe ich also folgendermaßen umgesetzt, indem ich überhaupt die Kraft garnicht erst berechne, sondern für jeden Körper die Beschleunigung. Die beiden Varianten sind zwei Seiten derselben Medaille, denn: berechne ich zunächst die für beide Objekte gültige Kraft, muss ich dazu beide Massen einbeziehen. Nachher, um die Beschleunigung des einzelnen Körpers zu errechnen muss ich aus der "gemeinsamen" Kraft die andere Masse wieder rausdividieren, denn in der Kraft stecken beide Massen und Beschleunigung ist
a = F/m (von F=m*a) Das hat mich dann zu dieser Variante geführt.
Deshalb ist auch meine "Accel" meiner Meinung auch genau das was sie bezeichnet nämlich die Beschleunigung und die wird bei mir in "DoStep" berechnet.


Der Funktion Move würde ich deshalb die Schrittweit als Parameter übergeben und den Code in der Funktion wie folgt anpassen:

freilich könnte man das machen aber genau wie vorher, muss man dann zweimal multiplizieren (mal 3 Komponenten) was ich mir durch den vorherigen Schritt erspart hab.


Man kann dann das Projekt als Ganzes kopieren und in der Methode Move die Midpoint-Methode implementieren.
Wie gesagt genau das wäre der nächste logische Schritt.
Diese Quelltextversion sollte nun genug optimiert sein um den Schritt zu gehen. Ausdrücklich ist jeder dazu eingeladen diese Version des Quelltext zu nehmen und Midpoint zu implementieren! :cool:


Eventuell sollte man die Methode Move auch umbenennen. Wie wäre es z.B. mit Planet::GetNextPosition(double schrittweite)?

einverstanden!


c) Bei der Klasse Planet fehlt mir noch der Destruktor, der den mit new reservierten Speicher für die Vektoren pos, Accel und Velocity wieder freigibt.
Gruß
Ja das ist richtig. War halt nicht so wild da die Planeten sowieso bis zum Ende vorhanden waren! Aber ist natürlich korrekt!

Andere Änderungen am Code:
- Planeten: Venus bis Saturn plus Erdmond
- Ausgabe geändert: Zeitschritte sind nun ca. ein Monat und nach 12 Schritten ( = 365 Tage) endet die Integration. Ausgabe nach jedem Monat. Ein Monat dauert bei mir nun 20sec.
- Klassen weiter aufgeräumt
- Genauigkeit: Fehler in X-Position nach 1 Jahr ca.102.000km (laut HORIZONS)
- Speedupdate aus letztem Beitrag von mir

Known Bugs:
- zu hohe Ungenauigkeit
- PSystem Destruktor fehlt (alle Planeten löschen) - unkritisch

Ab jetzt kommt in jedes Paket nurnoch die Quellcodedateien plus die geschwindigkeitsoptimierte EXE Datei.

RELEASE 02 (http://www.coding-facility.com/sr/PSystemC_02.rar)

fr0he 0stern

Bernhard
22.04.2011, 21:51
Das ist in dem Fall nicht korrekt da es sich um Zeiger handelt und auch wenn die Zeigervariable selbst gelöscht wird so bleibt der Inhalt des Speichers erhalten, weshalb ein delete erforderlich ist.
Hallo Stefan,

das delete sollte trotzdem raus aus der Schleife. Das frisst wirklich unnötig Rechenzeit. Man alloziert besser einmal außen den Speicher (am besten ein einziges Mal beim Programmstart) und beschreibt die vorhandenen Speicherbereiche immer wieder neu mit den temporären Werten. Beim Programmende wird der Speicher dann wieder freigegeben. Bei den paar double-Werten sollte man das gut hinbekommen.


Aus Gründen der "mathematischen Correctness" wäre das richtig, aber da bei dieser Variante die Gconst mit multipliziert werden muss, hab ich sie gleich beim Festlegen der Masse einbezogen. So spart man sich die eine Multiplikation dort wo es auf speed ankommt.
OK. Sehr gute Idee.


Ausdrücklich ist jeder dazu eingeladen diese Version des Quelltext zu nehmen und Midpoint zu implementieren! :cool:
Der Unterschied zwischen beiden Verfahren ist nicht sehr groß:

Euler:
x_{n+1} = x_n + h * v_n(x_n)
v_{n+1} = v_n + h * a(x_n)

Das ist momentan auch so implementiert. Die Addition wird dabei über die Methode PVektor::Add umgesetzt.

Bei der Mittelpunkts-Methode verwendet man dagegen:
x_{n+1} = x_n + h * v_n(x_n + h/2 * v_n(x_n))
v_{n+1} = v_n + h * a(x_n + h/2 * a(x_n))

Hast Du als Test mal nachgesehen in wieviel Tagen der Mond um die Erde kreist?
Gruß und frohe Ostern

Bernhard
22.04.2011, 22:20
Man alloziert besser einmal außen den Speicher (am besten ein einziges Mal beim Programmstart) und beschreibt die vorhandenen Speicherbereiche immer wieder neu mit den temporären Werten.
Ich sehe gerade, dass man dazu leider auch die Funktion PVektor::Sub umschreiben müsste, aber es lohnt sich. Du könntest z.B. die Adresse des Speichers für das Ergebnis als zusätzlichen Parameter übergeben und sparst dann pro Schleifendurchlauf das new und dann auch das delete.

SRMeister
22.04.2011, 22:49
Hallo Stefan,

das delete sollte trotzdem raus aus der Schleife.
Du hattest vollkommen recht. !!!!!! .
Die Zeit für 1 Monat ist von 21 auf 7 sec runter gegangen.
wow!


Hast Du als Test mal nachgesehen in wieviel Tagen der Mond um die Erde kreist?

nein, kommt bei Gelegenheit!

und Midpoint!

SR

Bernhard
22.04.2011, 23:19
Bei der Mittelpunkts-Methode verwendet man dagegen:
x_{n+1} = x_n + h * v_n(x_n + h/2 * v_n(x_n))
Hi Stefan,

mir fällt gerade noch ein, dass der Term v_n(x_n + h/2 * v_n(x_n)) dummerweise nicht direkt berechenbar ist. Der Algorithmus zur Mittelpunkts-Methode passt also wegen der zweiten Ableitung noch nicht vollständig.

Bitte also, gemäß meiner Anleitung, noch nicht implementieren.

Die Theorie muss erst passen, sonst ärgert man sich bloß, weil der Code nutzlos ist - Sch...
Gruß

Bernhard
23.04.2011, 02:04
Bei der Mittelpunkts-Methode verwendet man dagegen:
x_{n+1} = x_n + h * v_n(x_n + h/2 * v_n(x_n))
v_{n+1} = v_n + h * a(x_n + h/2 * a(x_n))
Ich habe da die Variablen x und v durcheinander gebracht. Richtig muss es heißen:

x_{n+1} = x_n + h * v_n + 0,5 * h^2 * a(x_n)
v_{n+1} = v_n + h * a(x_n + h/2 * v_n)

Die Beschleunigung muss bei diesem Verfahren pro Schritt, statt einmal, also zweimal berechnet werden und zwar an der Stelle x_n für den neuen Ort und an der Stelle x_n + h/2 * v_n für die neue Geschwindigkeit.

Wie man das herleitet kann man z.B. in den numerical recipes in C, Second edition, Kapitel 16.1 Runge-Kutta Method nachlesen.

SRMeister
23.04.2011, 12:32
Nur nochmal zum Verständniss,

x_n - Aktueller Ort
v_n - aktuelle Geschw.
a - Beschleunigung?
h - Schrittweite
n aktueller Zeitpunkt
n+1 nächster (zu berechnender) Zeitpunkt


die Beschl. muss also hier quasi zweimal nur als Zwischenprodukt berechnet werden.
deshalb muss ich wohl erstmal die berechnung derselben auch in Planet::SetNewPos() verlegen, wie du ja schonmal vorgeschlagen hast.
Ist der Teil hier
a(x_n + h/2 * v_n)
als mathematischer Term gemeint, oder ist damit die Beschleunigung a an der Stelle (x_n + h/2 * v_n) gemeint? Eher das Zweite sonst ergibts für mich wieder keinen Sinn.
Ansonsten aber schon.
Ich berechne a also mehrmals: einmal am Anfang also bei x_n und bei a(x_n + h/2 * v_n).



Euler:
x_{n+1} = x_n + h * v_n(x_n)
v_{n+1} = v_n + h * a(x_n)

Das ist momentan auch so implementiert. Die Addition wird dabei über die Methode PVektor::Add umgesetzt.

Da ist mir gerade was aufgefallen. Ich glaube bei mir ist das so nicht implementiert, sondern folgendermaßen:
v_{n+1} = v_n + a(x_n)
x_{n+1} = x_n + v{n+1}

1. fehlt das h in den Gleichungen, was aber durch TconsSquare = h*h bereits weiter vorne eingerechnet ist (möglicherweise ist frühe quadrieren falsch)
2. benutze ich v(n+1) um x(n+1) zu berechnen was möglicherweise auch falsch ist.

Ich denke zuerst müssen wir die Fehler da ausräumen, die möglicherweise für die Abweichungen verantwortlich sind.

Bernhard
24.04.2011, 11:32
x_n - Aktueller Ort
v_n - aktuelle Geschw.
a - Beschleunigung?
h - Schrittweite
n aktueller Zeitpunkt
n+1 nächster (zu berechnender) Zeitpunkt

"100 Punkte für den Kandidaten".

Bernhard
24.04.2011, 15:35
Ich habe da die Variablen x und v durcheinander gebracht. Richtig muss es heißen:

x_{n+1} = x_n + h * v_n + 0,5 * h^2 * a(x_n)
v_{n+1} = v_n + h * a(x_n + h/2 * v_n)

Die Beschleunigung muss bei diesem Verfahren pro Schritt, statt einmal, also zweimal berechnet werden und zwar an der Stelle x_n für den neuen Ort und an der Stelle x_n + h/2 * v_n für die neue Geschwindigkeit.
Hallo Stefan,

die Beschreibung der Formel a(x_n + h/2 * v_n) bedarf noch einiger Erklärung: Man muss bei der Berechnung der optimierten/zweiten Beschleunigung mit Hilfe der zwei Schleifen optimierte Planetenpositionen hernehmen. D.h. vorausgesetzt man kennt die Positionen und Geschwindigkeiten beim Schritt n, berechnet man daraus für jeden Planeten die optimierte Position r_i' = r_i + h/2 * v_i (wobei i die restlichen Planeten bezeichnet, die den Aktuellen beschleunigen) und benutzt diese N-1 Positionen für die Berechnung der Beschleunigung auf den aktuellen Planeten. Leichter geht es bei diesem Verfahren leider nicht.

Letztlich muss man also für jeden Zeitschritt sämtliche Positionen und Geschwindigkeiten für alle Planeten berechnen.
Gruß

Bernhard
24.04.2011, 18:42
Da ist mir gerade was aufgefallen. Ich glaube bei mir ist das so nicht implementiert, sondern folgendermaßen:
v_{n+1} = v_n + a(x_n)
x_{n+1} = x_n + v{n+1}

1. fehlt das h in den Gleichungen, was aber durch TconsSquare = h*h bereits weiter vorne eingerechnet ist (möglicherweise ist frühe quadrieren falsch)
2. benutze ich v(n+1) um x(n+1) zu berechnen was möglicherweise auch falsch ist.

Ich denke zuerst müssen wir die Fehler da ausräumen, die möglicherweise für die Abweichungen verantwortlich sind.
Hallo Stefan,

ich verwende aktuell die folgenden Änderungen für die Euler-Integration:

void Planet::GetNextPosition(double Tconst)
{
dx->x = Velocity->x;
dx->y = Velocity->y;
dx->z = Velocity->z;

// Neue Position
dx->Mul(Tconst);
pos->Add(dx);

// Neue Geschwindigkeit
Accel->Mul(Tconst);
Velocity->Add(Accel);

Accel->Reset();
}

und


void PSystem::DoStep()
{
steps++;
Planet *source;
Planet *target;
double dist;

for(source = sun; source != 0; source = source->next)
{
for(target = sun;target != 0;target = target->next)
{
if(source != target)
{
PVektor::Sub(tmp, target->pos, source->pos); // Beschleunigung berechnen
dist = tmp->Len();
tmp->Mul(target->mass/(dist*dist*dist));
source->Accel->Add(tmp); // Beschleunigung hinzufügen
}
}
}

for(source = sun; source != 0; source = source->next)
{
source->GetNextPosition( Tconst );
}
}

Bei der Ausgabe gebe ich die Geschwindigkeiten mit aus und bekomme dann nach einem Monat ganz geringfügig andere Werte bei der z-Position, bei v_y und v_z, aber die Integration müsste dann gemäß Theorie passen.

Bernhard
25.04.2011, 10:33
In der Funktion PSystem::AddObject habe ich ebenfalls die Multiplikation mit Tconst entfernt, damit in der Variable Velocity auch wirklich eine Geschwindigkeit gespeichert wird:


tmp->Velocity->x = vx;
tmp->Velocity->y = vy;
tmp->Velocity->z = vz;


Interessant wäre es noch Operatoren zu überladen, so dass man in der Klasse PVektor die Operatoren =, + und - benutzen kann, aber zuerst schreibe ich jetzt erst mal eine Version mit Midpoint. Du kannst mir gerne eine PN mit einer eMail-Adresse schicken. Ich schicke dann das Update mit dem neuen Code als zip-Datei.
MfG

SRMeister
26.04.2011, 02:45
Also dank Bernhard haben wir jetzt die Midpoint Methode.
Weiterhin sind jetzt alle Planeten incl. Erdmond in der Sim.
Der Fehler der ungenügenden Genauigkeit habe ich vermutlich gefunden: die Zeitfunktion die noch aus C-Zeiten stammt, rechnet anscheinend mit der Systemzeitzone (also unsere), jedenfalls beträgt der Fehler genau 6h.
Berücksichtigt man die Abweichung, beträgt der absolute Fehler der Position nach einem Jahr:
X = 730 km
Y = 5900 km
Z = 1640 km

Auch die Mondposition ist ähnlich gut.

Jedenfalls, wie gesagt, dank Bernhard haben wir nun die Midpoint Variante. Diese hat zunächst keine wirkliche Veränderung gebracht. Das heisst die Genauigkeit war vorher praktisch dieselbe, bezogen auf 1sec Schritte. Nun ist die Genauigkeit immernoch besser als vorher, wenn ich auf 100sec. Zeitschritte gehe. Da liegt der Vorteil dieser Methode.

Ich hoffe auf eine weitere Steigerung der Genauigkeit. Jedenfalls bringt die Einführung der Lichtgeschwindigkeit allein garnix, denn die Gravitation erfahren wir so als ob sie sich instantan ausbreitet.
Sicher gibt es irgendwelche linearen Approximationen für die ART Effekte, mal schauen.

Hier die aktuelle Version (http://www.coding-facility.com/sr/PSystemC_04.rar)

Stefan

Bernhard
26.04.2011, 05:44
In dem von euch vorgeschlagenen HORIZONS Katalog(Danke nochmal, genau was ich gesucht habe) werden perfekterweise die Objekte als Position und Geschwindigkeit angegeben.
Hier als Beispiel die Erde, Koordinatenursprung ist die Sonne, Einheit ist AU:
JDCT
X Y Z <-- Position
VX VY VZ <-- Geschwindigkeit
LT RG RR LT = Lichtzeit(day); RG = Distance from Center(AU); RR = Range Rate(Änderungsrate von RG pro Tag (AU/day)

2455652.500000000 = A.D. 2011-Apr-01 00:00:00.0000 (CT)
-9.855554845305842E-01 -1.864518242419004E-01 1.975847807252305E-05 <-- Position
2.933461796767275E-03 -1.697932543422671E-02 2.779216797989756E-07 <-- Geschwindigkeit
5.793060517472141E-03 1.003037335417702E+00 2.739049035701513E-04
Hallo Stefan,

mir ist momentan nicht klar, wie man an diese Zahlen kommt. Könntest Du bitte erklären, wie man das anstellt. Ich bekomme immer nur die Ephemeriden wie RA und Deklination in Abhängigkeit des Beobachtungsstandortes.

Die angegebenen Zahlen müssten sich aber auf den Schwerpunkt des Sonnensystems beziehen :confused:
MfG

Bernhard
26.04.2011, 05:51
Sicher gibt es irgendwelche linearen Approximationen für die ART Effekte, mal schauen.
Man könnte die Massen der Planeten trivialerweise über die relativistische Massenformel m_rel = m / sqrt(1-v^2/c^2) berechnen, da wir uns ja im Schwerpunktsystem befinden. Ich frage mich schon des längeren, ob man damit näherungsweise die Periheldrehung simulieren kann. Ich muss mir aber erst mal überlegen, wie sich das innerhalb der Midpoint-Methode implementieren läßt. Mit Runge-Kutta könnte man auch noch eine Version basteln, aber die relativistischen Effekte finde ich erst mal spannender.
MfG

SRMeister
26.04.2011, 11:13
Horizons Step by Step:

http://ssd.jpl.nasa.gov/horizons.cgi

Ephemeris Type : CHANGE ->
Vector Table

beim Rest sollte folgendes eingetragen sein:
Target Body [change] : Earth [Geocenter] [399]
Coordinate Origin [change] : Solar System Barycenter (SSB) [500@0]
Time Span [change] : Start=2012-03-31 05:00, Stop=2012-03-31 07:00, Step=10 m
Table Settings [change] : output units=KM-S; quantities code=2
Display/Output [change] : default (formatted HTML)

Leider bietet HORIZONS nur km als Einheit, was ich irgendwie "dumm" finde :) naja man muss also zum Exponenten immer 3 addieren :)
Ja und das mit der relativistischen Masse klingt interessant.

Bernhard
26.04.2011, 12:42
Wow, das ist ja wirklich äußerst komfortabel. Werde jetzt (auch) erst mal ein wenig "rumtesten". Danke.

Bernhard
26.04.2011, 18:57
Der Fehler der ungenügenden Genauigkeit habe ich vermutlich gefunden: die Zeitfunktion die noch aus C-Zeiten stammt, rechnet anscheinend mit der Systemzeitzone (also unsere), jedenfalls beträgt der Fehler genau 6h.
Das kann ich noch nicht nachvollziehen. Gerechnet wird mit 31536000 Sekunden = 365 Tage. Addiert man zum 01.04.2011 diese 365 Tage bekommt man - wie im Programm - den 31.03.2012.

Berücksichtigt man die Abweichung, beträgt der absolute Fehler der Position nach einem Jahr:
X = 730 km
Y = 5900 km
Z = 1640 km

Ich komme deswegen eher auf die folgenden Werte:

X = 109.459 km (Punkt als Trennzeichen für die Tausend)
Y = 640.785 km
Z = 30 km

Die Abweichung in Y ist größer als die Entfernung Erde - Mond und die Ergebnisse verbessern sich auch nicht, wenn man die Schrittweite von 10s auf z.B. 2s verkleinert.
Gruß

Bernhard
27.04.2011, 11:53
Man könnte die Massen der Planeten trivialerweise über die relativistische Massenformel m_rel = m / sqrt(1-v^2/c^2) berechnen, da wir uns ja im Schwerpunktsystem befinden. Ich frage mich schon des längeren, ob man damit näherungsweise die Periheldrehung simulieren kann.
Hallo Stefan,

angesichts des relativ großen verbleibenden Fehlers (den ich einfach mal auf Rechen- und Rundungsfehler zurückführe) fände ich einen Vergleich mit einer elliptischen Bahn interessanter. Man modelliert dabei die Bahn der Erde als Ellipse deren Brennpunkt im Schwerpunkt des Systems Erde-Sonne liegt.

Man braucht dann noch die Umlaufdauer, Bahnneigung der Ellipse, Länge des aufsteigenden Knotens und die Exzentrizität der Ellipse, um die Bahn eindeutig im dreidimensionalen Raum festzulegen. Man bekommt diese Daten auch mit Hilfe von HORIZONS mit der Einstellung

Ephemeris Type [change] : ELEMENTS

Natürlich sind diese Elemente aufgrund der Bahnstörungen durch die anderen Planeten nicht konstant, aber man könnte diese Werte einfach für einen Bahnumlauf als konstant annehmen und dann versuchen den resultierenden Fehler zu berechnen.

Um dabei auch beliebige Bahnabschnitte berechnen zu können, muss man noch wissen, wie schnell sich die Erde auf der Ellipse bewegt. Dafür gibt es aber eine Gleichung, die man dann ebenfalls numerisch lösen kann. Momentan weiß ich jedoch nicht, ob diese Gleichung bereits im Internet veröffentlicht wurde.
Gruß

Runzelrübe
27.04.2011, 13:45
Bitte versucht mal, die Anteile aufsteigend nach anzunehmender Größe zu addieren (in manchen Formeln durch simples Umstellen erreichbar).
Divisionen durch 2 sind zu ersetzen durch Multiplikationen mit 0,5.
Quadrate können durch Ausmultiplizieren einer eventuellen numerischen Approximation entgehen oder zumindest der Prüfung.
Wann immer eine Summe eine Multiplikation spart, anwenden (ja, ich habe im Beispiel auch den Pseudo-Code für a(...) umgestellt)

x_{n+1} = 0,5 * h * h * a(x_n) + h * v_n + x_n

oder

x_{n+1} = h * (0,5 * h * a(x_n) + v_n) + x_n

v_{n+1} = h * a(0,5 * h * v_n + x_n) + v_n

Die Änderungen der Positionen sind generell als kleiner zu erwartende Größen zu betrachten als die Position selbst. Ebenso Geschwindigkeitsänderungen. :)

---

G * Sum(G_free_part_i) wäre numerisch genauer als Sum(G * G_free_part_i), vor allem, wenn die Summe über die Anteile selbst auch wieder aufsteigend sortiert durchgeführt wurde. Also nicht die fertig ausmultiplizierten Anteile addieren sondern nach der Summenbildung noch eine Schleife über source und einmal die Multiplikation mit G für jede fertige Summe durchführen.

---

dist = tmp->Len() wird einzig und allein für dist*dist*dist benötigt.

Besser wäre, die lokale Variable dist fliegt ganz raus!

dist*dist*dist ersetzen durch tmp->LenSq()*tmp->Len()

oder durch

tmp->LenSq()*Sqrt(tmp->LenSq())

Somit wird die Abstandsberechnung auf dem Vektor nicht in ein lokales Hilfs-Double gepresst, nachdem sogar schon die Wurzel gezogen wurde. Selbst der Wurzelanteil fließt dann nur einmal ein.

LenSq() ist hierbei der quadratische Abstand vor dem Wurzelziehen (also bitte nicht die Wurzel ziehen und dann nochmal quadrieren ;)).
Len() kann für PVektor eventuell umgebaut werden als Sqrt(LenSq()).
Falls PVektor nicht verändert werden kann, eigene quadratische Längenberechnung einbauen oder ableiten.
Wäre eigentlich sogar toll, wenn hierbei die Summen der quadrierten Anteile nach Größe in LenSq einfließen (im Fall Solarsystem zumindest Y zuerst).

Die Objekte werden anhand ihrer Massen für das Array aufsteigend sortiert. Erst Asteroiden, dann Monde, dann Planeten, dann die Sonne.

Bernhard
27.04.2011, 15:55
Hi Rübe,

die Variable dist wird man sehr leicht los und mit der Schrittweite h kann ebenfalls sehr leicht etwas ökonomischer umgehen. Beide Punkte bringen zwar keinen Gewinn an Rechengenauigkeit (habe es gerade ausprobiert), machen aber den Code etwas übersichtlicher und schneller in der Ausführung. Danke also für denTip.

Wesentlich arbeitsintensiver ist der Umgang mit stark unterschiedlich großen Zahlen. Da kann man das Programm sicher noch stark verbessern. Gerade bei der Berechnung von Vektorlängen gibt es z.B. so Tricks wie: sqrt(x_1*x_1 + x_2*x_2) = fabs(x_1) * sqrt(1.0 + (x_2/x_1)*(x_2/x_1)), falls x_1 und x_2 ungefähr gleich groß sind.

Was das Ordnen von Summationsreihenfolgen bringen soll ist mir momentan nicht klar. Wenn bei einer Summe die größere Zahl um viele Zehnerpotenzen größer ist als die kleinere Zahl, wird die Kleinere eben irgendwann nicht mehr berücksichtigt und man hat einen Rechenfehler. Die Reihenfolge der Addition spielt da doch keine Rolle mehr, oder täusche ich mich da?
Gruß

Runzelrübe
27.04.2011, 16:45
Stell Dir vor, es würden maximal 4 Stellen in Deinen Datentyp passen und 8 könnten zusammengerechnet werden (interne Präzision).

Addiere nun mal 300 Zahlen, die so im Bereich 0,001 bis 999,9 liegen. Du kommst auf mathematisch 0,300 bis 299970,0.
4 Ziffern kannst Du nur im Typ abbilden. Daraus werden jetzt aber 0,300 bis 3,0E+5, denn 299970,0 wird auf 4 darstellbare Stellen gerundet und mit einem Exponenten versehen, dann in den Datentyp gepackt. Wir nehmen der Einfachheit halber mal die 999,9 für jede der 300 Zahlen.

Wurde die Summe jetzt einer Variable dieses Typs zugewiesen, haben wir die interne Präzision, die sich zumindest 299970,0 merken konnte, schonmal verloren.

Jetzt subtrahieren wir 4971,0 und addieren noch 1,4E+7. Einmal, bevor wir auf den Datentyp runden, einmal, nachdem wir das bereits gemacht haben. Dann kommen wir bereits auf 2 unterschiedliche Ergebnisse.

intern möglich: 299970,0 - 4999,0 + 14000000,0 = 14294971,0 -> 1,429E+7
bereits auf Datentyp: 3,00E+5 - 4999,0 + 14000000,0 = 14295001,0 -> 1,43E+7

Sieht ja nicht allzu wild aus. Die kleinen Zahlen hatten ja eh wenig Anteil. Das bisschen Rundung, naja, wer schreit danach. ;)

Okay, nun machen wir das mal in der "falschen" Reihenfolge.

Wir haben 1,4E+7 und addieren mal 300 Zahlen zwischen 0,001 und 999,9.
Wir machen das in einer Schleife, die diese Zahl immer wieder der Variablen zuweist (tmp = tmp + x)

Man mag nun denken, okay, kommt was ähnliches raus. Aber der Datentyp wird nach jeder internen Summe wieder zurückgecastet.

1,4E+7 + 0,001 -> 14000000,0 (nicht intern abbildbar, Größenordnung zu unterschiedlich [editiert, stand bis eben "noch abbildbar"]) -> 1,4E+7
1,4E+7 + 999,9 -> 14000999,9 -> 1,4E+7 (tja, wir haben eben nur 4 Stellen)

Nach 300 Additionen haben wir dann immernoch 1,4E+7.
Man könnte auch gern 1 Billiarde mal solche "kleinen" Zahlen zu 1,4E+7 addieren. Das Ergebnis würde und würde nicht größer werden, obwohl 1 Billiarde mal 999,9 insgesamt definitiv größer ist als 1,4E+7, nämlich 9,999E+14. Da sollte die 1,4E+7 irgendwann auf dem Weg zur finalen Summe ganz schön blöd aus der Wäsche schauen, aber sie ist immernoch da.

Nunja, ziehen wir mal jetzt noch die 4999,0 ab.

1,4E+7 - 4999,0 = 13995001 -> und wieder kommt 1,4E+7 raus (sowas aber auch ;))

Stellen wir mal 1,4E+7 der mathematisch genauen Summe von 14294971 (bei jeweils 999,9) gegenüber, haben wir einen Fehler von ca. 2,06%.
Bei der korrekten Reihenfolge und damit 1,43E+7 gegenüber 14294971 zumindest nur 0,03%.

Daher ist man auf der sichereren Seite, wenn man erst die kleinen Zahlen miteinander addiert und dann in der Größenordnung wächst. Oder aber, indem man, wie ich in einigen Posts zuvor erwähnte, Prädikate und Datentypen verwendet, denen die Größenordnungen egal sind und die ihre Fehleranteile kennen und weiterverwenden. Dann benötigt man auch keine Epsilonumgebungen dort, wo sie nicht hingehören und der Fokus kann auf dem Optimieren des Codes liegen. :)


-----------

Zum Thema Solarsystem: http://de.wikipedia.org/wiki/Sonnenmasse

Und die Einflüsse der Sonne sollten daher auch zuletzt betrachtet werden. Wir müssen im Computer manchmal erst die i-Tüpfelchen vor dem Offensichtlichen betrachten.

Gruß,
Rübe

Bernhard
27.04.2011, 17:31
Daher ist man auf der sichereren Seite, wenn man immer erst die kleinen Zahlen miteinander addiert und in der Größenordnung wächst. :)
So leuchtet es ein.
Danke.

SRMeister
27.04.2011, 20:17
Ich hab die Vorschläge von RR auch gerade durchprobiert.
Aber es ergibt sich tatsächlich keine wirkliche Verbesserung.

Meiner Meinung nach, ist die prinzipielle Art der Berechnung jetzt vollkommen korrekt, genau genug und nicht weiter zu verbessern. Abgesehen von relativistischen Effekten, die aber mit Sicherheit nicht für den Fehler verantwortlich sind!

Ich finde, dass die Position, wie ich bereits weiter oben sagte, nach einem Jahr um 6 Stunden falsch liegt. Es ist recht merkwürdig das es "genau" 6 Std. sind. Aber ich mache mittlerweile die Zeitfunktion nichtmehr verantwortlich dafür.

Ich vermute irgendeinen Fehler in den Vorgaben von Horizons.

Eine mögliche Erklärung wäre, dass sich das Sonnensystem selbst dreht gegenüber den Fixsternen, was aber nicht in den Koordinaten vorgesehen ist.

Ich habe auch durch die Hilfe von Horizons meine Vermutung weder bestätigen noch widerlegen können.
Aber:
Es gibt unter den Optionen bei Table Settings [Change] eine Option correction : -- selects level of correction to output vectors

dort ist momentan "none" eingegeben, es gibt aber 2 weitere Optionen:
- None ( Geometric States) (standart)
- Astrometric States
- Astrometric States corrected for stellar abberation

möglicherweise liegt hier der Fehler, denn "None" könnte falsch sein.

leider habe ich unter der Woche wenig Zeit um das zu testen. Ich komme aber drauf zurück falls mir niemand zuvorkommt ;)

Stefan

Runzelrübe
27.04.2011, 20:22
Ich hätt da noch einen für euch. Ich hab nämlich mein altes Programm ausgebuddelt, weil mich das Thema jetzt wieder angesteckt hat. :cool:

Anstatt sun -> pluto gegenüber sun -> pluto zu durchlaufen nehmen wir mal eine untere Dreiecksmatrix.

Bitte entschuldigt, wenn ich der Schnelle geschuldet unqualifizierten Code hier hinschreibe. Ich wollt das nicht vorher testen. Mein Programm läuft halt in C# und ich wollte euch zumindest den Ansatz vermitteln. Denn damit könnt ihr aus O(n^2) zumindest den kleinen Gauss rauskitzeln. :)

aus:




void PSystem::DoStep()
{
steps++;
Planet *source;
Planet *target;
double dist;

for(source = sun; source != 0; source = source->next)
{
for(target = sun;target != 0;target = target->next)
{
if(source != target)
{
PVektor::Sub(tmp, target->pos, source->pos); // Beschleunigung berechnen
dist = tmp->Len();
tmp->Mul(target->mass/(dist*dist*dist));
source->Accel->Add(tmp); // Beschleunigung hinzufügen
}
}
}

for(source = sun; source != 0; source = source->next)
{
source->GetNextPosition( Tconst );
}
}


mach:



void PSystem::DoStep()
{
steps++;
Planet *pI= sun;
Planet *pJ = sun->next;
double dist, dist3;
PVektor tmpI, tmpJ;
int iLen, i, j; //iLen mal bitte auf Anzahl Objekte setzen (wie gesagt, Code noch nicht angeschaut)!!!

for(i=0; i<iLen; i++)
{
for(j=i+1; j<iLen; j++)
{
PVektor::Sub(tmpI, pI->pos, pJ->pos);
PVektor::Sub(tmpJ, pI->pos, pJ->pos); //notwendig (doppelt gemoppelt), da Sub, Add usw. keine neue Instanz zurückgeben sondern die bestehende ändern (brauch aber zwei verschiedene Vektoren zum Ändern des jeweiligen Beschleunigungsanteils)
dist = tmpI->Len();
dist3 = dist*dist*dist;
tmpI->Mul(pJ->mass/dist3)
tmpJ->Mul(pI->mass/dist3)
pI->Accel->Sub(tmpI);
pJ->Accel->Add(tmpJ);

pJ=pJ->next;
}

pI=pI->next;
}

for(pI = sun; pI!= 0; pI=pI->next)
{
pI->GetNextPosition( Tconst );
}
}


Falls das jetzt nicht sofort jedem verständlich war: wir durchlaufen in der äußeren Schleife alle Objekte, in der inneren jedoch nur diejenigen, die noch vor uns liegen. Dafür rechnen wir aber die Auswirkungen beider Objekte gegeneinander in einem Schritt der innerern Schleife aus. Da derselbe Vektor entgegengesetzte Auswirkungen hat, subtrahieren wir in einem Fall und addieren im anderen. Tadaaaa.. anstatt bei 100 Objekten 100x100 Durchläufe mit insgesamt 10.000 Tests auf Gleichheit haben wir nur noch 5.050 Durchläufe ohne einen einzigen Test.

Viel wichtiger aber: wir haben immernoch das Problem mit der lustigen Sortierreihenfolge! Die Ungenauigkeiten bleiben leider erhalten. Lassen sich aber zumindest ein wenig wegoptimieren (hab ich eben nicht getan). :(

Runzelrübe
27.04.2011, 20:31
Es ist recht merkwürdig das es "genau" 6 Std. sind. Aber ich mache mittlerweile die Zeitfunktion nichtmehr verantwortlich dafür.
Ich vermute irgendeinen Fehler in den Vorgaben von Horizons.

Du machst Dir zuviele Sorgen.

Gehen wir mal die Sache ganz anders an. Was sind 6h. 1/4 Tag. Hmmmm...
Was ist alle 4 Jahre, aber nicht alle 100 Jahre, dafür aber alle 400 Jahre, doch wiederum nicht alle 11200 Jahre..... richtig :D

Und was ist 2012? Auch richtig! Nein, nicht die Maja, menno, denk ans Schaltjahr! :D

Ich nehme sehr stark an, dass der 29.02.2012 nur anteilig mit seinen 6h in die Positionsbestimmung eingeflossen ist. Probier mal bis 2013 zu rechnen und schau, ob dann 12h Unterschied sind. :) *Daumen drück*


[Ich werde nachher hoffentlich mal ein Video meiner alten Version posten - lasse sie derzeit für 2500 Sterne durchrechnen, 1s Film benötigt so ca. 30s Rechenzeit, ist leider (noch) nicht parallelisiert; die Ausgabe erfolgt per GDI unkomprimiert in ein AVI]

Bernhard
27.04.2011, 20:48
Ich nehme sehr stark an, dass der 29.02.2012 nur anteilig mit seinen 6h in die Positionsbestimmung eingeflossen ist.
Das is es nicht. Ich habe das mit einem Kalender nachgeprüft und der 29.02.2012 wird in der Rechnung zu 100% berücksichtigt.

Am liebsten würde ich alles mal mit 100 Stellen Genauigkeit rechnen und nicht mit double.

Runzelrübe
28.04.2011, 10:01
Am liebsten würde ich alles mal mit 100 Stellen Genauigkeit rechnen und nicht mit double.

Leider habe ich die Implementierung der Shewchuk (http://www.cs.berkeley.edu/~jrs/)'schen Prädikate (hier ein Link zur incircle-Implementierung in C (http://www.cs.cmu.edu/afs/cs/project/quake/public/code/predicates.c)) nur für .NET umgesetzt. Für Division, Wurzel, Potenzen, Trigonometrie muss man entweder sein Köpfchen noch selbst anstrengen oder aber gleich auf LEDA (http://www.mpi-inf.mpg.de/LEDA/) / CGAL (http://www.cgal.org/) zurückgreifen.
Viel Erfolg! Ich empfehle CGAL! Dort findet ihr auch viele andere Konzepte. :)

Bernhard
30.04.2011, 10:48
Ich empfehle CGAL!
Hallo Rübe,

vielen Dank für die interessanten Links, aber ich "verrate" einfach mal, dass ich momentan an einer eigenen Klasse bastle, die auf die Funktionen aus den "Numerical Recipes in C" zurückgreift. Mit diesen Funktionen habe ich bisher sehr gute Erfahrungen gemacht und in einer eigenen Klasse läßt sich damit ein sehr mächtiger Datentyp implementieren.

Mächtig deswegen, weil man in dieser Klasse das error handling selbst gestalten kann und dann auch Rundungsfehler getrennt analysieren kann. Der Exponent läßt sich in eine Variable vom Typ long auslagern. Das Vorzeichen kommt in eine Variable vom Typ short und bei der eigentlichen Zahl kann man dann angeben, wie viele Stellen vor und nach dem Komma berücksichtigt werden sollen.

Das Thema wird damit ziemlich komplex, aber immerhin sind Addition, Subtraktion und Multiplikation schon fertig. Für das N-Körper-Problem braucht man dann nur noch die Division und die Quadratwurzel, wobei beide Funktionen im Prinzip via recipes schon vorhanden sind. Man muss diese Funktionen dann nur noch in die vorhandene Klasse einbauen und Vorzeichen und Exponent getrennt verwalten.
MfG

SRMeister
01.05.2011, 20:39
Hi Leute :)

wow, endlich mal nen ganzen Tag Zeit gehabt und genutzt!

Erstmal die Fakten:

1. Habe Runge-Kutta-4 implementiert.
2. Fehler ist jetzt bei 1000sec Schrittweite ca. gleich wie bei Midpoint mit 100sec. Schrittweite. -> Extreme speed-up
3. Der Mathematische Fehler wird nicht besser, wir waren bereits mit Midpoint (100sec.) und Euler (1sec.) sehr nah dran.
4. Habe die Operatoren implementiert wie von Bernhard vorgeschlagen, außerdem alle PVektoren auf Variablen umgestellt (keine Zeiger mehr).

==> Der Fehler, mit dem wir hier zu kämpfen haben, ist nicht mathematischer Natur! Auch eine Verzehnfachung der Nachkommastellen wird daran nichts ändern!

==> Wir müssen andere physikalische Gesetzmäßigkeiten mit einbeziehen. Allg. RT. ODER der Fehler liegt in den fehlenden Massen (Jupitermonde... Oortsche Wolke...Kuipergürtel... Keine Ahnung) Gibt es Angaben dazu, wieviel Masse in dem Asteroidengürtel zwischen Erde, Mars und Jupiter existiert? Wenn da ein Erdäquivalent an Masse drin steckt, könnte das den Fehler schon erklären (Vermutung)

Ansonsten, nachwievor kommen wir an der Stelle raus, die Horizons am 31.03.2012 um 6:00 Uhr hat und nicht um 0:00 Uhr wie es korrekt wäre,. Das bedeutet, die Erde bewegt sich 0,0685% schneller als sie laut Horizons sollte.

Die Quellcodeversion basiert auf Bernhard seiner letzten,
allerdings ist sie massiv geändert wurden.
Euler und Midpoint funtkionieren nach wie vor, wurden aber ebenfalls an die gemachten Änderungen angepasst.
Wer jetzt noch Verfahren höherer Ordnung als RK4 einfügen will, braucht dazu fast keinen Aufwand mehr.

>> LINK (http://www.coding-facility.com/sr/PSystemC_06.rar) <<

Bernhard
02.05.2011, 08:57
1. Habe Runge-Kutta-4 implementiert.
Gratulation.

2. Fehler ist jetzt bei 1000sec Schrittweite ca. gleich wie bei Midpoint mit 100sec. Schrittweite. -> Extreme speed-up
Prima.

==> Der Fehler, mit dem wir hier zu kämpfen haben, ist nicht mathematischer Natur!
Die/Der erwähnte neue Klasse/Datentyp ist mittlerweile codiert, eingebaut und ganz grob getestet. Schwierigkeiten bereitet noch die Verwaltung des Exponenten, weil hier doch sehr oft sehr unterschiedlich große Zahlen miteinander verknüpft werden und Ziel ist es natürlich dabei möglichst wenig Informationen zu verlieren. Runzelrübe hatte das weiter oben bereits erwähnt. Ich werde der Sache also noch weiter nachgehen.

BTW: Ich habe neulich gelesen, dass es in PHP scheinbar ähnliche Funktionen gibt, die ebenfalls mit vorgebbarer Genauigkeit arbeiten.
Gruß

UMa
02.05.2011, 15:11
Hallo Stefan,

ich empfehle das von mir im Link von MAC (die weiteren Posts von mir in jenem Thread beachten) beschriebene Verfahren. Das sollte RK-4 um längen schlagen und liefert automische Fehlerabschätzung und Zeitschrittweitensteuerung. Zeitschritte um 1 Tag sollten für das Sonnensystem (mit Erdmond) möglich sein, ohne Genauigkeitsverlust.

Der Unterschied zwischen ART und Newton sollte für das Sonnensystem ca 1e-8 betragen,
Ich vermute daher noch einen Fehler. Evtl position von Erde und Mond bezüglich des Baryzenters? Oder die Massen sind nicht genau genug, man sollte immer dern Graviationsparameter GM nehmen (Genauigkeit 1e-11 oder so), nie(!) aus M und G ausrechen.

Ich würde empfehlen zunächt ein rein newtonsches N-Körperproblem zu implementieren.


Grüße

UMa

SRMeister
02.05.2011, 15:22
ja alles klar,
das hilft weiter.

Bis jetzt war das bestmögliche Ergebnis pro Verfahren bei allen Verfahren gleich, dh. das bestmögliche Endergebnis war bei jedem Verfahren immerwieder der selbe Zielort. Nur die Schrittgröße hat sich immerweiter erhöht bei vernachlässigbarem Fehler.

Ahja, das Programm setzt den GM selbst fest.
Habe aber bereits den GM Parameter bei Horizons entdeckt und werde den jetzt gleich mal testen!

Hoffen wir das Beste!

Zeitschritte um 1 Tag sollten für das Sonnensystem (mit Erdmond) möglich sein, ohne Genauigkeitsverlust.
Klingt auch gut! Schaue es mir mal an.

Stefan

Bernhard
02.05.2011, 15:24
ich empfehle das von mir im Link von MAC (die weiteren Posts von mir in jenem Thread beachten) beschriebene Verfahren.
Hallo UMa,

genau darum ging es doch die letzten paar Seiten? Ich dachte eigentlich das wäre die Midpoint-Methode??


Oder die Massen sind nicht genau genug, man sollte immer dern Graviationsparameter GM nehmen (Genauigkeit 1e-11 oder so), nie(!) aus M und G ausrechen.
hast Du dazu eine geeignete Quelle im Netz?
MfG

SRMeister
02.05.2011, 15:55
Hallo UMa und Bernhard,


hast Du dazu eine geeignete Quelle im Netz?

Die braucht er nicht, ich habs probiert, und DAS WARS !!!
Der Fehler ist sogut wie verschwunden !!!!

Ich glaube was UMa mit der anderen Methode meinte, war, den Fehler abzuschätzen, der ja mit abnehmender Schrittweite asymptotisch gegen Null läuft. Wenn man den Fehler bei großer Schrittweite kennt, braucht man ihn nur zum Ergebnis zu addieren, ohne wirklich die Schrittweite kleiner zu machen.
Es gibt auch einen Namen dafür, den hab ich schon irgendwo gelesen. Kann sein das es Leapfrog oder Verlet war.

Stefan

Nachtrag, der Fehler der Erde nach einem Jahr zu Horizons:
X = 9.131m
Y = 61.461m
Z = 345m

Betrag = 62,1 km oder 2 sekunden Bewegung

UMa
02.05.2011, 16:07
Hallo Stefan, hallo Bernhard,

Ich verwende "The JPL Planetary and Lunar Ephemerides, DE405/LE405"
http://iau-comm4.jpl.nasa.gov/relateds.html
da kann man sich eine Datei de405iom/de405iom.pdf herunterladen an deren Ende Massenverhältnisse und Positionen bzw Geschwindigkeiten angegeben werden.

Das ist nur das Basisverfahren.
http://www.astronews.com/forum/showthread.php?p=35493#post35493

Die entscheidenen Posts sind hier
http://www.astronews.com/forum/showpost.php?p=35612&postcount=145
Die Ergebnisse des Basisverfahrens werden extrapoliert. Das gesamte Verfahren ist auch als
http://en.wikipedia.org/wiki/Bulirsch-Stoer_algorithm
bekannt. Bei 7 Stufen d.h. 6facher Extropolation (d.h. Berechungen mit Schrittweiten von 1/1, 1/2, .. ,1/7) hat es Ordnung 14.
6 oder 7 Stufen sind für das Sonnensystem optimal.

Eine Ergänzung/Korrektur
http://www.astronews.com/forum/showpost.php?p=35729&postcount=147

Grüße

UMa

Bernhard
02.05.2011, 17:32
Die Ergebnisse des Basisverfahrens werden extrapoliert. Das gesamte Verfahren ist auch als
http://en.wikipedia.org/wiki/Bulirsch-Stoer_algorithm
bekannt. Bei 7 Stufen d.h. 6facher Extropolation (d.h. Berechungen mit Schrittweiten von 1/1, 1/2, .. ,1/7) hat es Ordnung 14.
6 oder 7 Stufen sind für das Sonnensystem optimal.

Eine Ergänzung/Korrektur
http://www.astronews.com/forum/showpost.php?p=35729&postcount=147
Vielen Dank UMa.

Genau diese Angaben hatten (mir) gefehlt und damit kann man dann das bestehende Programm sicher auch noch weiter verbessern. Supi :) .
MfG

Bernhard
02.05.2011, 17:48
Nachtrag, der Fehler der Erde nach einem Jahr zu Horizons:
X = 9.131m
Y = 61.461m
Z = 345m

Betrag = 62,1 km oder 2 sekunden Bewegung
Super. So langsam werden die Werte interessant!
MfG

Bernhard
05.05.2011, 16:25
Hi Stefan,

ich will jetzt mal nicht zuviel versprechen, aber es sieht so aus, als wäre die Klasse BigFloat fertig. Ich habe eben eine Rechnung mit 60 Nachkommastellen gestartet, um zu sehen wie lange eine Simulation über einen Monat mit der Midpoint-Methode dauert und komme auf eine Rechenzeit von knapp 10 Stunden. Es wäre also durchaus interessant den Bulirsch-Stoer-Algorithmus zu programmieren, um die Rechenzeit nochmal so stark zu senken, dass man dann auch die erhöhte Genauigkeit nutzen kann.

Möchtest Du noch weitere Versionen veröffentlichen oder "stinkt" Dich das Thema mittlerweile eher an?? Für mich alleine würde ich das Projekt jetzt erst mal auf CD brennen und archivieren. Wenn Dich UMas Algo aber noch interessiert, schreibe ich auch noch gerne Folgeversionen. Ob ich dabei aber die neue Klasse ebenfalls als Freeware herausgebe, weiß ich noch nicht. Unter Umständen ist das aber auch gar nicht nötig.
Gruß

SRMeister
05.05.2011, 18:38
Hi Bernhard,
nein nein, mich stinkt das garnicht an!
Das Bulirsch Stoer Verfahren möchte ich zu gerne noch implementieren.
Ebenso auch relativistische Korrekturterme.
Ich bin halt Selbstständig und hab unter der Woche immer wenig Zeit zum Programmieren. Sowas mach ich nur gerne "am Stück" also wenn ich mir nen ganzen Tag zeitnehmen kann.

Ich, und ich vermute die Anderen Interessenten hier, freuen uns doch über jede gemachte Verbesserung. Ich hoffe wir machen hier zusammen weiter!

Ich habe da noch eine Idee: Wenn du statt des Midpoint Verfahrens, die RK4 Variante nimmst, kannst du bei gleicher Genauigkeit etwa 10mal größere Zeitschritte gehen (die allerdings jeweils 4 mal länger dauern, aber netto ergibt das trotzdem einen Zeitgewinn).
Ebenso würde ich vorschlagen, wir gestalten den Bulirsch Stoer Code so, dass man zwischen Midpoint u. RK4 einfach wechseln kann.
Willst du den machen? Ansonsten mach ich den vllt. am WE

Grüße

Bernhard
05.05.2011, 21:02
Ebenso auch relativistische Korrekturterme.
Die oben geäußerte Idee mit der relativistischen Planetenmasse habe ich vorhin noch mit double-Genauigkeit und Midpoint-Methode gerechnet. Die Korrekturen sind aber leider minimal und bringen IMO das Programm nicht wirklich weiter. Die Änderungen am Code kann ich Dir bei Bedarf aber gerne zuschicken.

Ich hoffe wir machen hier zusammen weiter!
Gern.

Ebenso würde ich vorschlagen, wir gestalten den Bulirsch Stoer Code so, dass man zwischen Midpoint u. RK4 einfach wechseln kann.
Willst du den machen? Ansonsten mach ich den vllt. am WE
Da bist Du dann voraussichtlich schneller als ich. Ich habe heute erst mal den Wikipedia-Artikel dazu gelesen, weil ich dieses Verfahren bisher noch nie verwendet habe. Aktuell könnte ich nur die hier ( http://de.wikipedia.org/wiki/Richardson-Extrapolation ) beschriebene Richardson-Extrapolation einbauen.
Gruß

Bernhard
07.05.2011, 00:33
Hallo Stefan,

ich bin gerade mit der Durchsicht der neuesten Version (6) fertig. RK4 ist sehr übersichtlich programmiert -> Super.

Die Funktion PSystem::Evaluate läßt sich noch etwas beschleunigen:


if( i == 0 )
{
for(source = sun; source != 0; source = source->next)
{
source->der[0].dv = 0;
for(target = sun; target != 0; target = target->next)
{
if(source != target)
{
tmp = target->Pos - source->Pos;
tmp *= target->mass / (tmp.Len()*tmp.SqrLen());
source->der[0].dv += tmp;
}
}
source->der[0].dx = source->Vel;
}
}
else
{
for(source = sun; source != 0; source = source->next)
{
source->der[i].dv = 0;
for(target = sun; target != 0; target = target->next)
{
if(source != target)
{
tmp = (target->Pos + target->der[i-1].dx * dt) - ( source->Pos + source->der[i-1].dx * dt) ;
tmp *= target->mass / (tmp.Len()*tmp.SqrLen());
source->der[i].dv += tmp;
}
}
source->der[i].dx = source->Vel + source->der[i-1].dv * dt;
}
}

Die Fallunterscheidung wird dann pro Aufruf nur noch einmal statt N(N-1) ausgeführt. Bei i=0 kann man auch die unveränderte Differenz der Positionen nehmen und spart damit noch mal einige Rechenoperationen.

Bei der Berechnung der Länge der Vektoren würde ich eventuell wieder die einfache Funktion sqrt(x^2+y^2+z^2) verwenden. Die von mir vorgeschlagene Version stammt noch aus Zeiten, wo der Exponent auf 38 begrenzt war. Da macht dann das Quadrieren bei großen Zahlen Probleme, die man bei double und unseren Zahlen aber nicht mehr hat.
MfG

Bernhard
07.05.2011, 10:37
Hallo Stefan,

bevor Du die nächste Version baust, würde ich auch noch die Funktionen Planet::SetTrueMass und Planet::GetTrueMass rausnehmen. Die Gravitationskonstante kann ebenfalls komplett raus und die Variable für die Planetenmasse kann man dann z.B. GM nennen.

Bei den Startwerten muss man dann natürlich statt der Masse GM (aus Horizons) einsetzen:


sys->AddObject("Sonne",1.32712440018E+20,
-6.080091705285974E+08, 2.454333708238292E+07, 2.193321432139981E+06,
2.809726412820165E+00, -1.036343441959919E+01, -4.204860196381065E-02 ); // Sonne

sys->AddObject("Merkur",22032.09E9,
-5.331602529851193E+10, 1.424323914756870E+10, 6.000210156549231E+9,
-2.277058581046753E+04, -4.495485791176564E+04, -1.582633849867808E+03); // Merkur

/*
sys->AddObject("Ceres",63e9,
3.450197517424105E+11, -2.719974973036755E+11, -7.217640948809467E+10,
1.023158320287534E+04, 1.294765961349763E+04, -1.480644559407215E+03); // Ceres (1)

sys->AddObject("Juno",1882148849.4,
-4.217785409007722E+11, 3.611805220272993E+10, 8.817606494305907E+09,
-5.567212093261015E+03, -1.582354566781498E+04, 3.813944367844741E+03); // Juno (3)

sys->AddObject("Pallas",14.3e9,
1.288986305799774E+11, -4.038017005052425E+11, 2.682923164674586E+11,
1.422726738202094E+04, 1.371255413474907E+03, -2.139616981343534E+03); // Pallas (2)

sys->AddObject("Vesta",17.8e9,
1.756133416483675E+10, -3.228910857346381E+11, 7.507208978128448E+09,
2.092921321574529E+04, 5.360629186823918E+02, -2.559040856144347E+03); // Vesta (4)
*/

sys->AddObject("Venus",324858.63E9,
2.397967758924109E+10, -1.059717457755222E+11, -2.868972017005123E+09,
3.388294521475002E+04, 7.783094129350907E+03, -1.848596675496016E+03); // Venus

sys->AddObject("Erde",398600.44E9,
-1.474370019336122E+11, -2.789279589304088E+10, 2.955826247744262E+06,
5.079162483215040E+03, -2.939896910566963E+04, 4.812093925554706E-01); // Erde

sys->AddObject("Mond",4902.798E9,
-1.470514277349460E+11, -2.801509364690417E+10, 3.730129039033502E+07,
5.382156089456466E+03, -2.847733344864282E+04, 2.235354455693894E+01); // Mond

sys->AddObject("Mars",42828.3E9,
2.035953497078640E+11, -3.456857769669849E+10, -5.736961786047079E+09,
4.976094385716297E+03, 2.594977520891352E+04, 4.218104801061298E+02); // Mars

sys->AddObject("Jupiter",126686511E9,
7.115603707882032E+11, 2.013825854433948E+11, -1.677091861778326E+10,
-3.715952626588193E+03, 1.319564076159422E+04, 2.830947269851247E+01); // Jupiter

sys->AddObject("Saturn",37931207.8E9,
-1.396843425883871E+12, -3.381973875783179E+11, 6.146826298318255E+10,
1.753564063629373E+03, -9.410322369209542E+03, 9.449519324100386E+01); // Saturn

sys->AddObject("Uranus",5793966E9,
3.003955040895567E+12, 2.614531534364399E10, -3.882080718055668E+10,
-1.089777236983593E2, 6.492185531504123E+03, 2.552972796165465E1); // Uranus

sys->AddObject("Neptune",6835107E9,
3.826621024919128E+12, -2.346728069113412E+12, -3.986079005702972E+10,
2.805108135465661E+03, 4.664771972182217E+03, -1.611737087776297E2); // Neptune


Wie Du siehst, habe ich auch mal die vier größten Asteroiden berücksichtigt. Die können aber den verbleibenden Fehler von rund 60 km auch nicht erklären.

Interessanterweise liegen die Rundungsfehler bei der Midpoint-Methode ungefähr in dieser Größenordnung. Die Rechnung mit einer Genauigkeit von 60 Stellen über knapp 12 Tage (also vom 01.04.2011 bis 12.04.2011) lieferte bereits eine Korrektur von gut 3 km.

Ich werde als nächstes also mal versuchen Runge-Kutte mit 60 Stellen zu rechnen.
Beste Grüße

EDIT: So richtig toll wäre noch die Umrechnung der kartesischen Koordinaten in RA und Deklination. Dann könnte man die berechneten Positionen auch mit Beobachtungsdaten vergleichen. In diesem Zusammenhang verstehe ich noch nicht, woran die Achsen des kartesischen Koordinatensystems festgemacht werden. Der Ursprung ist gleich dem Schwerpunkt, das ist klar. Aber in welche Richtung zeigen eigentlich die Achsen? Nachdem die Erde laut Horizons auch in z-Richtung eine Bewegung hat, liegt die x- und y-Achse nur näherungsweise in der Ekliptik.

Bernhard
11.05.2011, 22:11
Ich habe da noch eine Idee: Wenn du statt des Midpoint Verfahrens, die RK4 Variante nimmst, kannst du bei gleicher Genauigkeit etwa 10mal größere Zeitschritte gehen (die allerdings jeweils 4 mal länger dauern, aber netto ergibt das trotzdem einen Zeitgewinn).
Ebenso würde ich vorschlagen, wir gestalten den Bulirsch Stoer Code so, dass man zwischen Midpoint u. RK4 einfach wechseln kann.
Willst du den machen? Ansonsten mach ich den vllt. am WE
Hallo Stefan,

ich hoffe Du bist wohlauf, weil es wieder ein paar Neuigkeiten zu dem Programm gibt. Ich habe RK4 mal für einen Monat mit 60 Stellen Genauigkeit gerechnet, dabei aber nur unwesentliche Korrekturen < 1 km gefunden. Verringert man den Zeitschritt bei RK4 von 1000s auf 100s bringt das auch keine deutliche Verbesserung, was bedeutet, dass uns das Bulirsch-Stoer-Verfahren ebenfalls erst mal nicht weiter bringt.

Daran anschließend fragt man sich natürlich woher denn nun die verbleibenden 60 km Fehler kommen. Ein Blick in das Buch "Allgemeine Relativitätstheorie" von T. Fließbach hilft aber etwas weiter, weil dort der Wert für die Periheldrehung der Erde mit 5'' pro Jahrhundert angegeben ist und das sind gut 36 km pro Jahr. Allerdings dreht diese Periheldrehung in die "falsche" Richtung und kann deswegen den Unterschied zu den Werten vom NASA-Server nicht unmittelbar erklären. Trotzdem zeigt die Größenordnung dieser Korrektur IMO, dass wir so langsam an die Grenzen des ursprünglichen Konzeptes stoßen.
MfG

Runzelrübe
12.05.2011, 10:12
Bitte versucht mal, eure Zeitschritte nicht auf 100s oder 1000s zu setzen sondern rechnet nur mit echten Zweierpotenzen wie 128 oder 1024. Dadurch muss die Mantisse der Gleitkommazahl nicht verändert werden und nur der Exponent erhöht sich. Abgesehen vom Geschwindigkeitsvorteil verursacht ihr so weniger Fehler auf den ULPs, unabhängig von der Genauigkeit.

SRMeister
12.05.2011, 23:15
Die oben geäußerte Idee mit der relativistischen Planetenmasse habe ich vorhin noch mit double-Genauigkeit und Midpoint-Methode gerechnet. Die Korrekturen sind aber leider minimal und bringen IMO das Programm nicht wirklich weiter. Die Änderungen am Code kann ich Dir bei Bedarf aber gerne zuschicken.

Interessant ist das allemal. Würde gerne einen Codeschnipsel dazu lesen!



Da bist Du dann voraussichtlich schneller als ich.
Wie man sieht, sieht man nicht viel von mir im Moment :(
Es gibt im Moment keine Freizeit für mich...

Folgende Vorschläge von dir:


Die Funktion PSystem::Evaluate läßt sich noch etwas beschleunigen:

Bei der Berechnung der Länge der Vektoren würde ich eventuell wieder die einfache Funktion sqrt(x^2+y^2+z^2) verwenden.

würde ich auch noch die Funktionen Planet::SetTrueMass und Planet::GetTrueMass rausnehmen. Die Gravitationskonstante kann ebenfalls komplett raus und die Variable für die Planetenmasse kann man dann z.B. GM nennen.


Wie Du siehst, habe ich auch mal die vier größten Asteroiden berücksichtigt.

Finde ich insgesamt sehr interessant, ich beauftrage dich hiermit die Version 07 zu machen, die du ja sicher eh schon hast.
Schick sie mir per EMail und ich lad sie auf den Server.



Die Rechnung mit einer Genauigkeit von 60 Stellen über knapp 12 Tage (also vom 01.04.2011 bis 12.04.2011) lieferte bereits eine Korrektur von gut 3 km.
3km sind immerhin 5%. Wenn es in die richtige Richtung läuft, wär das nicht übel. Aber ich denke mit der Rechengenauigkeit brauchen wir uns nicht vor Horizons zu verstecken, viel genauer haben die es sicher auch nicht gemacht.
Leider findet man auch in der Horizons Anleitung kaum Infos zur Genauigkeit.
Der Fehler liegt sicher wiedermal woanders.



EDIT: So richtig toll wäre noch die Umrechnung der kartesischen Koordinaten in RA und Deklination. ... In diesem Zusammenhang verstehe ich noch nicht, woran die Achsen des kartesischen Koordinatensystems festgemacht werden.
Das steht in der Horizons Doku.
Die XY- Ebene ist die Ebene, die die Erde und die Sonne aufspannten und zwar (ich glaube) am 1.1.2000. Die nennen das Epoche J2000 und das ist sozusagen ein Referenzzeitpunkt.
Die positive Z Richtung wurde einfach definiert: sie ist die Normale auf der XY-Ebene, die auf der Seite liegt, auf der auch der Erd-Nordpol liegt. Etwas umständlich formuliert, aber der Erd Nordpol zeigt eben nicht in genau dieselbe Richtung wie +Z(ist nicht parallel).



Ich habe RK4 mal für einen Monat mit 60 Stellen Genauigkeit gerechnet, dabei aber nur unwesentliche Korrekturen < 1 km gefunden. Verringert man den Zeitschritt bei RK4 von 1000s auf 100s bringt das auch keine deutliche Verbesserung, was bedeutet, dass uns das Bulirsch-Stoer-Verfahren ebenfalls erst mal nicht weiter bringt.
Naja, immerhin können wir die Zeitschritte immer größer machen. Das bringt uns, dass wir langfristige Untersuchungen machen können, bspw. wie stabil sind Planetenbahnen, bspw. bei kleinen Modifikationen oder der Einfluss dieses hypothetischen Doppelsterns der Sonne.
Da spielt dann auch die Genauigkeit nicht mehr eine so große Rolle, statistisch sollten die Fehler sich ja aufheben(auf lange Sicht gesehen)


Daran anschließend fragt man sich natürlich woher denn nun die verbleibenden 60 km Fehler kommen.
Ich habe da schon eine Vermutung:
Sicher genau wie ich, hast du die GM-Werte aus Horizons genommen. Allerdings hat UMa etwas weiter vorn einen Link gepostet, der auch die GM Werte angibt, aber in diesem Dokument stehen andere Werte. Obwohl Horizons angeblich diese Werte zugrunde gelegt hat.... Also irgendwas stimmt da nicht. Vielleicht gibt Horizons aktuellere oder ältere Werte auf ihren Seiten an, mit denen sie aber garnicht gerechnet haben.

http://iau-comm4.jpl.nasa.gov/relateds.html
... Datei de405iom/de405iom.pdf


Bitte versucht mal, eure Zeitschritte nicht auf 100s oder 1000s zu setzen sondern rechnet nur mit echten Zweierpotenzen wie 128 oder 1024.
Noch ein kleines ToDo für Bernhard und die 07er Version ;)

Ich lese immer mit nur hab ich im Moment selten die Zeit die ich mir für ordentliche Beiträge nehmen möchte. Ich find es aber echt großartig dass es weitergeht!
Außerdem, wo bleiben die bisher immer hilfreichen Ideen der stillen Mitleser?

Stefan

Bernhard
13.05.2011, 00:54
Interessant ist das allemal. Würde gerne einen Codeschnipsel dazu lesen!
Version 07 enthält die Methoden PSystem::Evaluate_Relativistic und PSystem:: DoStepRungeKutta_Relativistic. Dort wird ganz naiv die relativistische Planetenmasse eingesetzt. Die Korrekturen, die sich daraus ergeben liegen im Bereich 5 Meter.

Schick sie mir per EMail und ich lad sie auf den Server.
Done.

SRMeister
13.05.2011, 10:36
Jup, hier (http://www.coding-facility.com/sr/Code_07.rar)haben wir die aktuelle Version.

Runzelrübe
13.05.2011, 13:28
Ich habe mich eben mal ein wenig mit der 3D Darstellung dieser Szenerie beschäftigt und hätte da mal einen Wunsch. Könntet ihr eure Ergebnisse jeder Berechnung eventuell in eine CSV Datei gießen (id, x, y, z reichen aus je Objekt), so dass ich eure Ergebnisse (und auch Verbesserung/Verschlechterung selbiger) mal visuell darstellen kann? Danke!

Da ich ein Pathtracing eingebaut habe, würde sogar einer visuellen Gegenüberstellung zur Horizons Datanbank nichts im Wege stehen. Manchmal erkennt das menschliche Auge Fehler anhand von Bewegungen bei übertrieben hohen Geschwindigkeiten. Wenn die Daten also nicht erst live berechnet werden müssen, kann die Darstellung auch viel zügiger ablaufen. Systematisch größer werdende Abweichungen wären nicht mehr nur eine Zahl.

-----

Während der Betrachtung der Pfade ist mir eine Idee gekommen. Mal angenommen, wir würden mit einem Teleskop Langzeitbelichtungen derjenigen Sternsysteme durchführen, in denen wir bereits bestätigte Transits haben. Dann könnten wir mit diesen Aufnahmen einen weit vom Zentralgestirn entfernten Planeten immer dann erkennen, wenn dieser aus unserer Sicht seine horizontale Richtung wechselt. Die Wahrscheinlichkeit sollte ja eigentlich hoch sein, dass die Ebene durch den bestätigten Transit auch für andere Planeten gilt.

Das sähe dann ungefähr so aus:

S... Stern
P... Planetenumkehr


PPP.............................................S. ........................................PPP

Bei PPP sollte während der Langzeitbelichtung ein durchaus markanter Punkt zu sehen sein. Hintergründe könnte man ja anhand von Kurzzeitbelichtungen subtrahieren.

Wäre das machbar? Sinnvoll? Wird das eventuell schon gemacht?

sirius3100
13.05.2011, 14:27
Ich verstehe nicht genau was gemeint ist: Wer soll seine Richtung ändern? Der Transitplanet, oder der Stern, oder der Planet selbst?
Und warum zwangsläufig horizontal und nicht auch vertikal (radial schließe ich mal aus)?
Am ehesten hört sich das ganze ja nach Astrometrie an, aber dafür brauch ich natürlich keinen Transitplaneten in dem System. Das "funktioniert" auch ohne. Astrometrie hat halt den Nachteil das sie explizit von unserem Abstand zum beobachteten Stern abhängt.
Oder soll der Planet an den markierten Punkten der Bahn direkt beobachtet werden?
Direkte Beobachtung ist natürlich umso leichter je heller das Objekt leuchtet und je weiter es vom Zentralgestirn weg ist. Funktioniert aber deshalb eigentlich nur für schwere Gasriesen (andere Objekte sind einfach nicht hell genug um sie direkt zu sehen), und auch bei denen muss man tricksen.

Bewegt
13.05.2011, 14:54
S... Stern
P... Planetenumkehr


PPP.............................................S. ........................................PPP

Bei PPP sollte während der Langzeitbelichtung ein durchaus markanter Punkt zu sehen sein. Hintergründe könnte man ja anhand von Kurzzeitbelichtungen subtrahieren.

Wäre das machbar? Sinnvoll? Wird das eventuell schon gemacht?
geniale Idee
die funktioniert um so besser, je mehr sich der Beobachter auf Rotationsebene befindet

Runzelrübe
13.05.2011, 15:46
Wer soll seine Richtung ändern? Der Transitplanet, oder der Stern, oder der Planet selbst?

Ein (noch nicht entdeckter) Planet innerhalb eines Sternsystems, für das bereits ein Transitplanet bekannt ist.


Dann könnten wir mit diesen Aufnahmen einen weit vom Zentralgestirn entfernten Planeten immer dann erkennen, wenn dieser aus unserer Sicht seine horizontale Richtung wechselt.


Und warum zwangsläufig horizontal und nicht auch vertikal (radial schließe ich mal aus)?

Das habe ich sehr unglücklich ausgedrückt. In Systemen, in denen es aus unserer Sicht heraus Transitplaneten gibt (wir haben nachgewiesen, wie sie vor ihrem Zentralgestirn vorüberziehen), schauen wir aller Wahrscheinlichkeit nach nahezu senkrecht auf den Normalenvektor der Rotationsebene. :)

Dadurch dürften wir eine Planetenbahn als Linie betrachten (und bitte keine Spitzfindigkeiten jetzt!), auf denen der Planet von unserem Standpunkt aus oszilliert. Den Rest der Erklärung kannst Du jetzt gern nochmal durchlesen.

sirius3100
13.05.2011, 16:03
Hört sich also doch nach Direct Imaging an. Bisher konnten halt nur Exoplaneten mit einem Mindestabstand von ca 8AU (Beta Pictoris b) von ihrem Zentralgestirn direkt abgebildet werden.

Dein Vorschlag läuft also darauf hinaus bei der Auswahl zu beobachtender Sterne solche vorzuziehen bei denen bereits ein Transitplanet entdeckt wurde, da sich damit die Chance erhöht dass ein Planet von uns aus betrachtet weit genug vom Zentralgestirn weg ist um ihn direkt beobachten zu können.

Habe ich das jetzt richtig verstanden? Oder ist ganz was anderes gemeint?

Runzelrübe
13.05.2011, 16:26
Hast Du richtig verstanden. Deine Einwände bezüglich Planetengröße und Albedo sind natürlich nicht ohne - daher auch meine Frage, ob dieser Gedanke denn so abwegig sei. Da ich sehe, dass sich daraus eine eigenständige Diskussion entwickelt, die den ursprünglichen Inhalt des Themas überdeckt, habe ich das Thema "Weitere Planeten in Systemen mit Transitplaneten erkennen (http://www.astronews.com/forum/showthread.php?t=5194)" aufgemacht. :)

UMa
13.05.2011, 16:44
Hallo Stefan, hallo Bernhard,

ich habe mal mein altes Programm ausgegraben.
Ich habe vom Startpunkt de405iom.pdf Datei genau 360 Tage sowohl relativistisch als auch mit Newton vorrausberechnet. Der Unterschied im Ort der Erde ist 49791 m. Mit Horizons habe ich noch nicht vergleichen können.
Des weiteren könnte es Unterschiede im Zeitablauf geben. Horizons nimmt vermutlich Erdzeit, während ich in der Zeit eines weit entfernten Beobachters rechne (denke ich). Vielleicht erklärt das zusammen die ca. 60km Abweichung.

Nur relativistische Massen oder so einzusetzen verbessert das nicht, ihr müsstet eine bessere Approximation an die ART rechen als Newton z.B. die Einstein-Infeld-Hoffman Gleichungen.
http://en.wikipedia.org/wiki/Einstein%E2%80%93Infeld%E2%80%93Hoffmann_equations
Ergänzung: Die dort verlinkte erste Referenz
Standish, Williams. Orbital Ephemerides of the Sun, Moon, and Planets
düfte eigentlich alles weitere erklären.

Die GM Werte aus der Datei sind vielleicht für die Planeten mit Monden, die aus Horizons ohne Monde?!

Ich habe feste Zeitschritte von 1, 2, 3, 4, 5 und 6 Tagen verwendet. Die Unterschiede sind kleiner als zwischen relativistisch als auch mit Newton, zwischen 1d und 2d weniger als 1 m in der Erdposition.
Die Rechenzeiten für 1000 Jahre (genauer 360000 Tage) liegen bei wenigen Sekunden bis einer Minute je nach Version und Zeitschrittweite. Der Unterschied zwischen 1d und 2d Zeitschrittweite nimmt dabei für der Erde auf über 6km nach ca. 1000Jahren zu.

Ich würde nicht erwarten das RK4 + Extrapolation besser ist als Bulirsch Stoer, könnt ihr aber natürlich selbst ausprobieren.

Grüße

UMa

Bernhard
14.05.2011, 00:45
Ich würde nicht erwarten das RK4 + Extrapolation besser ist als Bulirsch Stoer, könnt ihr aber natürlich selbst ausprobieren.
Hallo UMa, Hallo Stefan,

ich habe jetzt doch mal das Bulirsch-Stoer-Verfahren gemäß UMas Anleitung eingebaut. Man bekommt bei damit bei achtfacher Schrittweite und gleicher Rechenzeit praktisch die gleichen Ergebnisse, wie mit RK4.

BTW: In den numerical recipes wird die Extrapolation etwas anders mit Hilfe von Polynomen, bzw. rationalen Funktionen durchgeführt. UMas Anleitung ist da wesentlich einfacher zu programmieren.

Die EIH-Gleichungen sind noch in Bearbeitung.
MfG

Bernhard
14.05.2011, 14:41
Nur relativistische Massen oder so einzusetzen verbessert das nicht, ihr müsstet eine bessere Approximation an die ART rechen als Newton z.B. die Einstein-Infeld-Hoffman Gleichungen.
http://en.wikipedia.org/wiki/Einstein%E2%80%93Infeld%E2%80%93Hoffmann_equations
Hallo UMa,

der Tip ist wirklich Spitze. Mit diesen Gleichungen kann man die Abweichung von etwas mehr als 60 km auf etwa 3 km reduzieren und dieser Wert liegt dann größenordungsmäßig in der Größe der Zeitdilatation zwischen Erdzeit und der Zeit im Schwerpunktsystem :cool: .


Ergänzung: Die dort verlinkte erste Referenz
Standish, Williams. Orbital Ephemerides of the Sun, Moon, and Planets
düfte eigentlich alles weitere erklären.
Sehr interessant, dass dort die Bewegung der Sonne aus der numerischen Integrationsroutine herausgenommen wird.


Die GM Werte aus der Datei sind vielleicht für die Planeten mit Monden, die aus Horizons ohne Monde?!
Genau so ist es scheinbar.
Beste Grüße

Runzelrübe
14.05.2011, 17:06
@SRMeister:

Habe mir eben Deinen Code angeschaut und noch ein paar geschwindigkeitsfördernde Optimierungsvorschläge, die Dir zumindest bei Deiner Umsetzung der relativistischen Berechnung helfen werden.

Aus:



PSystem::PSystem(double Timestep)
{
c = 299792458.0;
Tconst = Timestep;
sun = 0;
}
....
void PSystem::AddObject(char *Name, double Mass, double x, double y, double z, double vx, double vy, double vz)
{
...


mach:



PSystem::PSystem(double Timestep)
{
c = 299792458.0;
invC2 = 1.0/(c*c);
Tconst = Timestep;
tConstHalf = 0.5*Tconst;
sun = 0;
}
....
void PSystem::AddObject(char *Name, double Mass, double x, double y, double z, double vx, double vy, double vz)
{
numObjects +=1;


Dann kannst Du aus:



/c/c


das hier machen:



*invC2


und aus:



Evaluate(Tconst * 0.5,


das hier:



Evaluate(tConstHalf ,


Einer der wirklichen Geschwindigkeitsbooster, den ich Dir immernoch ans Herz lege, ist die Halbierung Der benötigten Schleifen je Zeitschritt:

Aus:



for(source = sun; source != 0; source = source->next)
{
source->der[0].dv = 0;
for(target = sun->next; target != 0; target = target->next)
{
if(source != target)
{


mach:



int i,j;
source = sun;

for(i=0; i<numObjects; i++)
{
source->der[0].dv = 0;
source = source ->next;
}

source = sun;
target = sun->next;

int outerLoopCount = numObjects-1;
for(i=0; i<outerLoopCount; i++)
{
target = source ->next;
for(j=i+1; j<numObjects; j++)
{
... Anteil wie bisher auf source anwenden
... umgekehrter Anteil (andere Masse und Richtung) auf target anwenden
target = target->next;
}
source = source->next;
}



Du hast nur die halbe Anzahl Durchläufe, sparst Dir die n^2 Prüfungen, ob source != target und musst auch nicht den Vektor, den quadratischen Abstand, den Abstand (Wurzelberechnung!) und den inversen Abstand (Division!) doppelt berechnen. Alles, was dazu kommt sind vier Multiplikationen und drei Vorzeichenwechsel (dazu wird ein Bit umgesetzt, ist ein Witz).

Bernhard
14.05.2011, 20:51
Du hast nur die halbe Anzahl Durchläufe
Hallo Rübe,

ich habe Stefan gerade die neueste Version mit sämtlichen Tips von UMa geschickt. Die Variablen invC2 und tHalfConst habe ich gleich miteingebaut. Die Anzahl der Durchläufe läßt sich IMO nicht reduzieren, weil sämtliche Körper in der Simulation unterschiedliche Massen haben. Man hat damit n^2 verschiedene Beschleunigungen, die immer wieder komplett neu berechnet werden müssen. Eine Umkehrung der Richtung der Beschleunigung alleine reicht dabei nicht aus.

Eventuell hat UMa seine Software noch besser als die neue Version 08 implementiert, weil die scheinbar um einiges schneller läuft.

Nach Einfügen der Körper Pluto, Ceres, Pallas, Vesta und Juno benötigt die Bulirsch-Stoer-Integration der EIH-Gleichungen bei einer Schrittweite von 8640 Sekunden und einer Integration über knapp fünf Jahre, d.h. bis zum 30.03.2016 00:00:00, etwa zweieinhalb Minuten Rechenzeit auf einem aktuellen Standard-PC. Der absolute Fehler in der Position der Erde beträgt im Vergleich zu den Werten aus HORIZONS etwa 53 km.
MfG

Bernhard
14.05.2011, 22:26
Der absolute Fehler in der Position der Erde beträgt im Vergleich zu den Werten aus HORIZONS etwa 53 km.

Berücksichtigt man auch noch den gravitativen Zeitdilatationseffekt der Erde, kann man den Fehler auf ca. 51 km reduzieren :) .

Runzelrübe
14.05.2011, 22:46
Die Anzahl der Durchläufe läßt sich IMO nicht reduzieren, weil sämtliche Körper in der Simulation unterschiedliche Massen haben. Man hat damit n^2 verschiedene Beschleunigungen, die immer wieder komplett neu berechnet werden müssen. Eine Umkehrung der Richtung der Beschleunigung alleine reicht dabei nicht aus.

Hi Bernhard,

geht wohl. Schau Dir mal folgenden funktionierenden C#-Code an (Real ist hierbei der nach Shewchuk implementierte Datentyp - dynamische Präzisionserhaltung garantiert):

Stell doch einfach mal den Code zweier Körper gegenüber. Dann siehst Du, dass Du nur (G)Masse austauschen und den Vektor umdrehen musst.. et voila.



public void CalculateAccelleration()
{

Real xx,yy,zz,sqrDistance,distance,invDistance3;
Real multiplierIII,multiplierJJJ;
IMovableMass objIII, objJJJ;

var objects = this.Formation.AsEnumerable().ToArray();
var objectsCount = objects.Length;
for (int i = 0; i < objectsCount; i++)
{
objects[i].AX = 0;
objects[i].AY = 0;
objects[i].AZ = 0;
}

var objectsCount1Less = objectsCount - 1;
for (int i = 0; i < objectsCount1Less; i++)
{
objIII = objects[i];
for (int j = i + 1; j < objectsCount; j++)
{
objJJJ = objects[j];
xx = objJJJ.X - objIII.X;
yy = objJJJ.Y - objIII.Y;
zz = objJJJ.Z - objIII.Z;
sqrDistance = xx * xx + yy * yy + zz * zz;
distance = Real.Sqrt(sqrDistance);
invDistance3 = 1.0 / (sqrDistance * distance);
multiplierIII = objJJJ.Mass * invDistance3;
multiplierJJJ = objIII.Mass * invDistance3;
objIII.AX += multiplierI * xx;
objIII.AY += multiplierI * yy;
objIII.AZ += multiplierI * zz;
objJJJ.AX -= multiplierJ * xx;
objJJJ.AY -= multiplierJ * yy;
objJJJ.AZ -= multiplierJ * zz;
}
}
}

Bernhard
15.05.2011, 08:20
geht wohl.
Hallo Rübe,

Danke für den Code. Ich werde mal versuchen diesen Trick auch auf die Integration der EIH-Gleichung anzuwenden. Die Integration dort läuft aktuell schon etwas langsam.
MfG

Runzelrübe
15.05.2011, 13:21
Ich werde mal versuchen diesen Trick auch auf die Integration der EIH-Gleichung anzuwenden. Die Integration dort läuft aktuell schon etwas langsam.

Hi Bernhard,

viel Erfolg und auch Spaß dabei! :)

Ich hab da noch ein paar Vorschläge. Lass Dich aber bitte nicht gleich davon verwirren. Für den Anfang reicht auch erstmal die Reduktion der Schleifendurchläufe insgesamt, ganz ohne zusätzliche Vorbereitungen.

Ich würde Dir sonst nämlich gern ein mehrstufiges Herangehen empfehlen. Erst die Komponenten errechnen, dann zusammenfügen.

Häufig verwendete Faktoren würde ich mir in lokale Variablen packen.

invC2
invSevenHalfC2

Nach dem Ansatz untere Dreiecksmatrix kannst Du zuerst einmal die benötigten Einzelkomponenten ausrechnen. Ich würde vorschlagen, ein struct dafür einzubauen und dieses je gerichteter Körperkombination in ein lokales Array zu speichern. Wichtig ist, dass diess Array beim Zugriff nach AB und BA unterscheiden kann. Die Indizierung ist Dir selbst überlassen.

Die sofort berechenbaren Einzelkomponenten sind:

nAB, nBA
vA2, vB2
vAB, vBA
rAB, r2AB, invRAB, invR2AB

Ich würde noch zwei weitere Variablen einführen:

pAB = xA-xB (x soll hier jeweils der Ortsvektor sein)
invR3AB = 1.0/(r2AB*rAB)

Somit können die Stellen, an denen "nBA/r2AB" steht, in "pBA*invR3AB" umgewandelt werden und die Stellen, an denen "xA-xB" steht, in "pAB".

Jetzt könnte man hingehen und auch noch die Terme in Variablen packen. Wenn Du das tust, dann lass Dich noch nicht von der Formel irritieren, die nach A, B und C unterscheidet. Du benötigst in den Vorbereitungen die Kombinationen AB und BA. Mehr nicht. Denn auch AC und BC werden in der Schleife irgendwann durchgenommen und aufgesammelt. Erst durch den Zugriff beim Zusammenrechnen greifst Du auf die Logik zurück, jeweils das Körperpaar A, B mit C bis N zu verbinden.

Sieht nach verdammt viel Aufwand aus, aber Deine Zugriffe beim Zusammenrechnen können das um Längen wieder reinholen.

Beim Kombinieren haben wir ja drei ineinander verschachtelte Schleifen.

for i in [0: (N-2)]
for j in [(i+1): (N-1)]
for k in [(j+1): N]

Wenn wir an alle voneinander abhängigen Terme per Zugriff rankommen, schicken wir die Berechnungsliste ja nur noch in die Pipeline. Somit wird dieser Punkt auf einmal der schnellste. Falls Dir das jetzt nicht sofort etwas sagt, hier nochmal ein kleines Beispiel:

Die drei Zeilen:

a = 3*x+1
b = 4*y+1
c = 5*z+1

können gleichzeitig ausgerechnet werden, sogar parallel (wenn man es integriert).

Jedoch:

a = 3*x+1
b = 4*y+a <- Zack, Unterbrechnung, wir benötigen a, um weiterzumachen.


So, da bin ich jetzt tief genug hinabgetaucht. Wie Du es im Endeffekt schneller machst, welche Variablen Du nimmst und wieviel Aufwand Du reinstecken möchtest für das schnellste Ergebnis, das ist Dir überlassen. Und Schritt für Schritt ist da immernoch die beste Lösung. :)

Bernhard
16.05.2011, 10:03
Hi Rübe,

die nichtrelativistische Rechnung wird durch die Berücksichtigung von \vec{r}_ab = -\vec{r}_ba schon mal um einen Faktor 2 (!) schneller und ist damit ein wirklich wichtiger Hinweis. Die relativistische Rechnung profitiert davon dann auch, weil die Beschleunigungen auf der rechten Seite der EIH-Gleichungen einfach aus der nichtrelativistischen Rechnung genommen werden.

Tricks, wie die Einführung von Konstanten wie SevenHalfC2 bringen ebenfalls Rechenzeit und machen den Code auch lesbarer. Bleibt dann noch die etwas schwierigere Strukturierung der Additionen, damit man auch hier die gleiche Arbeit nicht mehrfach macht.

Man kann mit der aktuellen und noch nicht öffentlich zugänglichen Version mit Runge-Kutta locker mal über 10 Jahre rechnen und das bei einer Schrittweite von 100s. Speziell bei dieser Rechnung waren die Fehler im Vergleich zur Rechnung über ein Jahr wieder deutlich kleiner (?). Auf jeden Fall macht es ziemlich Spaß mit dem Programm und HORIZONS etwas herumzuspielen.
Beste Grüße

SRMeister
16.05.2011, 12:05
Hier ist die aktuelle Version 08! (http://www.coding-facility.com/sr/Code_08.zip)
btw. good work, keep going

SR

Nachtrag:

und schon gehts hier (http://www.coding-facility.com/sr/Code_09.zip)weiter!

Bernhard
16.05.2011, 17:23
keep going
Hallo SR und R.Rübe,

Version 09 ist schnell genug, um 3650 (knapp zehn Jahre) Tage mit einer Schrittweite von 8640 Sekunden zu rechnen. Die nichtrelativistische Rechnung nach Bulirsch-Stoer wird bei mir in ca. 12s gerechnet. Die relativistische Rechnung braucht etwa 94s, ist aber leider von den Ergebnissen her schlechter, als die nichtrelativistische Rechnung.

Wer das Programm testen will: In der Datei PSystem.h kann in Zeile 7 ein Define gesetzt werden, das über einige Direktiven (#ifdef, ...) den Code steuert. Mit diesem Define wird relativistisch gerechnet. Wird die Zeile dagegen auskommentiert (und das Programm neu) übersetzt, wird nichtrelativistisch gerechnet. Im Code an SR ist die Zeile auskommentiert.
MfG

Bernhard
16.05.2011, 20:44
ist aber leider von den Ergebnissen her schlechter, als die nichtrelativistische Rechnung.
Bei einer Rechnung über 20 Jahre ist die relativistische Rechnung wieder deutlich besser.

Es zeigt sich dabei auch, dass die Windows-Kalenderfunktion noch falsch angesteuert wird. Erstens sollten bei den verwendeten Parametern (Vielfache von 86400) keine Uhrzeiten, wie 23:00 entstehen und zweitens kann man das korrekte Datum mittels des julianischen Datums (s. HORIZONS-Ephemeriden) nachprüfen. Bei der Rechnung über zwanzig Jahre beträgt dabei der Fehler vom ersten Hinsehen einen ganzen Tag.

Vorsicht also mit den Kalenderangaben beim Ausprobieren.
MfG

SRMeister
17.05.2011, 02:13
So ich habs endlich geschafft, eine Windows Version fertigzustellen.
Die zeitkritischen Funktionen (ich sagmal alles was innerhalb von DoSteps() passiert) werden in einer C++ DLL ausgeführt und sind somit so schnell wie die C++ .EXE datei. Die Windowsoberfläche ist C# basiert. Der Wrapperteill (Die Schnittstelle zwischen C++ und .NET) ist in C++.NET geschrieben. Es sind somit 3 Programmiersprachen enthalten.
Das ganze unterteilt sich dann logischerweise in mehrere Projekte:
PSystemDLL <-- enthält alle C++ bestandteile, basierend auf der Version 09 von Bernhard. Es wurde lediglich eine CountPlanets() Fkt. hinzugefügt da sie notwendig wurde.
PSystemWrapper Ist in C++.NET geschrieben und importiert die PSystenDLL in .NET
PSystemWrapperTest Ist eigentlich nicht nötig, aber war zum Testen gedacht. Es ist eine Konsolenanwendung, geschrieben in C# und stellt damit das DotNetGegenstück zur alten PSystemC.exe dar
NewSimulator ist die Windows Anwendung. Sie ist eine .exe und benötigt die PsystemDLL.dll und PSystemWrapper.dll im gleichen Ordner zum starten. Sie besteht wiederum aus einem Hauptprogramm und einem Benutzersteuerelemt "SimViewControl" was die Darstellung übernimmt und somit auch einfach in andere Programme eingebunden werden kann.

VisualStudio 2010 ist für den Quellcode benötigt.

Hier gehts zur 32-bit Binary (http://www.coding-facility.com/sr/WinSim09_Binary.rar)
Hier zum Quellcode (http://www.coding-facility.com/sr/WinSim09_Source.rar)

Have Fun!

Bernhard
17.05.2011, 20:07
Hi SR,

die grafische Oberfläche ist ein wirklich nettes Feature und die kleinen Kringel mit und ohne Farbe machen sich IMO ebenfalls recht gut. Trotzdem hätte ich noch ein paar Wünsche:

a) So weit ich weiß, wird unser Planetensystem üblicherweise von oben gezeigt, so dass die Planeten gegen den Uhrzeigersinn, also in mathematisch positivem Sinn, drehen.

b) Die vier größten Asteroiden Ceres, Pallas, usw. gehören von der Position her eigentlich zwischen Mars und Jupiter. Diesen Fehler habe ich in das Programm gebracht, weil es bei der Konsolenanwendung eigentlich keine Rolle spielt. Bei der grafischen Darstellung ist die korrekte Reihenfolge dagegen schon wichtig, weil sich der Normalbenutzer lieber geordnet von innen nach außen klicken will.

c) Ist es eventuell möglich mit der Programmmaximierung (zweiter Button ganz rechts oben, links neben dem Schließen Button: "x") auch das Display zu vergrößern. Man müsste dann nicht mehr mit der Maus scrollen, sondern könnte die Größe des Displays auch über die Größe des Programmfensters einstellen. Die Bedienbuttons unten in Programm sollten dabei auch einen festen Abstand zum Fensterrand haben.

Ich denke, diese drei Punkte würden die Bedienung des Programmes eventuell noch etwas komfortabler machen. Ansonsten finde ich die Oberfläche sehr gut gelungen.
MfG

Bernhard
18.05.2011, 15:25
Horizons nimmt vermutlich Erdzeit
Hallo UMa,

Horizons rechnet mit CT, was laut Manual für "coordinate time" steht. Vermutlich wird also bei Horizons, genau wie bei Deinem Programm, die Zeit im Schwerpunktsystem (= weit entfernter Beobachter) ohne weitere Korrekturen verwendet.

Unschön finde ich, dass diese Zeit bei Horizons nur indirekt (??) mit der UT zusammenpasst. Die Erdgeschwindigkeit durchschreitet beispielsweise den Nullpunkt am 1.1.2000 nicht um 11:58:xx sondern um 10:10:xx.
MfG

Bernhard
21.05.2011, 09:23
Hallo Stefan,

ich habe inzwischen eine neue Version für die Berechnung des N-Körper-Problems geschrieben und dabei alles etwas besser aufgeräumt. Sämtlicher Speicher, der auf dem Heap temporär angelegt wird, wird jetzt mit dem Programmende auch wieder freigegeben. So wird das von guter Software generell erwartet. Version 09 gibt die verlinkte Liste nicht korrekt frei. Der Speicherplatz für die einzelnen Körper wird deswegen nicht korrekt freigegeben.

Aus der Klasse PSystem habe ich die drei Klassen CNewton, CNewtonSP und CEIH gemacht. CNewton berechnet einfach das N-Koerper-Problem mit den bisher beschriebenen Methoden Euler bis Bulirsch-Stoer a la UMa. Bei der Klasse CNewtonSP wird die Bewegung der Sonne getrennt berechnet und zwar so, dass der Schwerpunkt des gesamten Systems konstant bleibt. Der wackelt nämlich bei Klasse CNewton ziemlich herum. Bei einer Rechnung über ein Jahr wandert der z.B. rund 35 km.

Bei einer Rechnung über 10 Jahre ist die Rechnung mit CNewtonSP leider ungenauer, als die Rechnung mit CNewton, was zeigt, dass auch die neue Version für Präzisionsrechnungen noch immer nicht taugt. Zwischen Bulirsch-Stoer und Runge-Kutta gibt es im Ergebnis keine nennenswerten Unterschiede.

Hoffe die neue Version gefällt (s. eMail). Kommentare dazu sind stets willkommen.
MfG

Bernhard
23.05.2011, 20:50
Inzwischen habe ich noch zwei Programmierfehler in der Berechnung mit erhöhter Genauigkeit (60 Nachkommastellen) gefunden. So konnte ich die Simulation über einen Monat mit Runge-Kutta bei einer Schrittweite von 1000 Sekunden sowohl mit double und der erhöhten Genauigkeit rechnen.

Der Unterschied liegt dabei interessanterweise im mm-Bereich. Im beschriebenen Fall sind die Rundungsfehler also zu vernachlässigen. Von den 15 Stellen der double-Ergebnisse sind also 13 korrekt.
MfG

UMa
25.05.2011, 09:22
Hallo Bernhard, hallo Stefan,

habt ihr mal ausprobiert wie sich der Fehler bei den einzelnen Verfahren als Funktion der Schrittweite verhält? Ihr könntet einen doppelt logarithmischen Plot für den Fehler (der Erdposition) in Abhängigkeit von der Zeitschrittweite machen, erstmal für die newtonschen Gleichungen. Als Referenz erstmal eine Rechnung mit kurzer Zeitschrittweite wählen.
Dann könntet ihr feststellen, welche Zeitschrittweite ihr noch wählen könnt ohne das der Fehler zu groß wird. 1000 Sekunden kommt mir noch etwas kurz vor, ich hatte ja min meinen Verfahren über einen Tag wählen können.

Grüße

UMa

Bernhard
25.05.2011, 13:30
Hallo UMa,

Stefan meldet sich seit einigen Tagen nicht mehr und ich vermute, dass er gerade im Urlaub verweilt oder auch eine gewisse Auszeit nimmt.

Mir persönlich erscheinen momentan auch die Werte von Horizons aus verschiedenen Gründen als ungeeignet und möchte deswegen als nächstes eine Umrechnungsfunktion in RA und Deklination einbauen, damit man die berechneten Werte mit Programmen wie Cartes Du Ciel vergleichen kann.

Generell würde ich Ephemeridenrechnung sowieso eher anders betreiben, z.B. als Zweikörperproblem mit Störtermen. Es gibt dazu, wie gesagt, das schöne Buch von O. Montenbruck im Verlag Sterne und Weltraum. Für längere Simulationszeiträume sind solche Verfahren einfach wesentlich übersichtlicher und netter zu programmieren.
MfG

Bernhard
25.05.2011, 19:35
1000 Sekunden kommt mir noch etwas kurz vor, ich hatte ja min meinen Verfahren über einen Tag wählen können.
Hi UMa,

die 1000 Sekunden hatte ich, so wie Stefan, eher für Runge-Kutta verwendet. Bei Bulirsch-Stoer bekommt man mit 8640 Sekunden (1/10 Tag) praktisch die gleichen Werte, wie bei Runge-Kutta mit 1000 Sekunden. Dass Bulirsch-Stoer ziemlich unabhängig von der verwendeten Schrittweite ist, hatte ich bei verschiedenen, kleineren Tests auch gesehen. Der Code sollte also bis hierher halbwegs passen :) .
MfG

SRMeister
26.05.2011, 21:37
So habe die grafische Version mal endlich weiter verbessern können. Finde sie macht jetzt eindeutig mehr her, habe die genannten Umschläge umgesetzt plus was mir noch so eingefallen ist. Wem fällt mehr ein?

In der Binär-Version sind 3 kompilierte Dateien, plus eine C++ Datei. Die C++ Datei muss zu Bernhards Quellcode hinzugefügt werden, dann muss das ganze Projekt einfach als DLL Datei kompiliert werden und kann dann die entsprechende DLL Datei der Binärdateien ersetzen. So kann jeder ohne Kenntnisse des Windows- Teils, also mit reinen grundlegenden C++ Kenntnissen theoretisch den "Core"-Quellcode weiter verbessern. Um DLL Sachen usw. braucht man sich also nicht kümmern.

Bei Bedarf lade ich auch gern wieder den kompletten Quellcode des Windows Projektes auf den Server.

Binary_10 (http://www.coding-facility.com/sr/Binary_10.rar)

Code_10_2 (http://www.coding-facility.com/sr/Code_10_2.rar) wie mir Bernhard gesendet hat, plus die oben genannte .cpp datei um es als .dll zu kompilieren

grüße,
Stefan

SRMeister
27.05.2011, 01:03
Version 09 gibt die verlinkte Liste nicht korrekt frei. Der Speicherplatz für die einzelnen Körper wird deswegen nicht korrekt freigegeben.
Das ist super und werden wir bald brauchen, wenn man die einzelnen Verfahren gegeneinander antreten lässt und den Fehler direkt auswerten lässt.



Bei der Klasse CNewtonSP wird die Bewegung der Sonne getrennt berechnet und zwar so, dass der Schwerpunkt des gesamten Systems konstant bleibt. Der wackelt nämlich bei Klasse CNewton ziemlich herum. Bei einer Rechnung über ein Jahr wandert der z.B. rund 35 km.
Ich denke das ist aber kein Fehler sondern einfach der RV-Effekt oder? Du siehst ja den Schwerpunkt nicht in den Daten, sondern nur den Nullpunkt. Oder bleibt der Schwerpunkt tatsächlich konstant?

Überhaupt, ich könnte spaßeshalber RV Kurven ausgeben lassen, fällt mir dabei ein!



Bei einer Rechnung über 10 Jahre ist die Rechnung mit CNewtonSP leider ungenauer, als die Rechnung mit CNewton, was zeigt, dass auch die neue Version für Präzisionsrechnungen noch immer nicht taugt. Zwischen Bulirsch-Stoer und Runge-Kutta gibt es im Ergebnis keine nennenswerten Unterschiede.
Was war denn die genaueste Methode bisher? Du hattest doch mal geschrieben, du hattest den Fehler auf unter 10km nach 1Jahr ?


Hoffe die neue Version gefällt (s. eMail). Kommentare dazu sind stets willkommen.
Muss unbedingt noch mehr testen, was die Verfahren angeht! :(



Stefan meldet sich seit einigen Tagen nicht mehr und ich vermute, dass er gerade im Urlaub verweilt oder auch eine gewisse Auszeit nimmt.
Ja das war korrekt. Es gibt da im Moment einige Probleme. Die Auszeit ist zum größten Teil leider unfreiwillig. Dann ist es schwer, sich komplett neu reinzuarbeiten. Allein das RK umzusetzen hat 2 ganze Tage hohe Konzentration gekostet, du kennst das ja sicher. Wenn man dann nur hin und wieder mal ne Stunde findet, fängt man auch lieber nicht mit solchen Monstern an.


Mir persönlich erscheinen momentan auch die Werte von Horizons aus verschiedenen Gründen als ungeeignet und möchte deswegen als nächstes eine Umrechnungsfunktion in RA und Deklination einbauen, damit man die berechneten Werte mit Programmen wie Cartes Du Ciel vergleichen kann.

Schau mal in die .exe rein. Ich denke viel fehlt da nichtmehr, um zu RA und Dec. von der Erde aus, zu kommen.

Generell würde ich Ephemeridenrechnung sowieso eher anders betreiben, z.B. als Zweikörperproblem mit Störtermen.
Also die Bahnberechnung per Kepler, plus Störterme. Soll sehr gut funktionieren, wie du sagst und nicht nur übersichtlicher sondern auch schneller zu sein. Aber ich denke das kommt dann wieder nicht an die von uns bzw. dir erreichte Genauigkeit ran. Ist mehr was für CUDA und 10.000 Körper.

Bei Bulirsch-Stoer bekommt man mit 8640 Sekunden (1/10 Tag) praktisch die gleichen Werte, wie bei Runge-Kutta mit 1000 Sekunden.
Wir müssen unbedingt mehr über die Genauigkeiten rausfinden, wie UMa auch sagt. Die wichtigste Frage ist hier, wir müssen uns ein Kriterium für eine Mindestgenauigkeit finden, und dann rausfinden, welches Verfahren diese Genauigkeit "am schnellsten" erreicht, also nicht abhängig von der Schrittgröße, sondern von der CPU Zeit. Es muss da irgendwo ein Optimum geben, was wir finden müssen. (Für eine vorgegebene Genauigkeit)

Da lässt sich in C# auch einiges viel leichter erreichen als auf C++ Konsole, denn du hast da die ganzen Windows Controls, wie Buttons, Tabellen und Charting für die Fehlerplots, was einfach die Bedienung weiter erleichtert.

Stefan

Bernhard
28.05.2011, 01:02
Das ist super und werden wir bald brauchen, wenn man die einzelnen Verfahren gegeneinander antreten lässt und den Fehler direkt auswerten lässt.
Hi Stefan,

schön, dass Du Dich meldest! Die Aufteilung in verschiedene Klassen ergab sich natürlich genau aus dem Wunsch die Verfahren über möglichst kleine Code-Änderungen gegeneinander antreten zu lassen.


Ich denke das ist aber kein Fehler sondern einfach der RV-Effekt oder?
Wofür steht RV?


Du siehst ja den Schwerpunkt nicht in den Daten, sondern nur den Nullpunkt. Oder bleibt der Schwerpunkt tatsächlich konstant?
Der Schwerpunkt ist zu 1.0 / Gesamtmasse * Summe (\vec{r}_Planet * m_Planet) definiert. Somit kann innerhalb der Klasse CNewton nach jedem Zeitschritt der Schwerpunkt des gesamten Systems mit einer zusätzlichen Schleife über alle Körper berechnet werden. Für den 01.04.2011 bekam ich über dieses Verfahren eine Abweichung des Schwerpunktes vom Koordinatenursprung von etwa 35 km.


Was war denn die genaueste Methode bisher? Du hattest doch mal geschrieben, du hattest den Fehler auf unter 10km nach 1Jahr ? Ich habe bei der neuesten Version Runge-Kutta mit CNewton als Sieger in Erinnerung, müsste das aber nochmal konkret nachsehen.


Wir müssen unbedingt mehr über die Genauigkeiten rausfinden, wie UMa auch sagt. Die wichtigste Frage ist hier, wir müssen uns ein Kriterium für eine Mindestgenauigkeit finden, und dann rausfinden, welches Verfahren diese Genauigkeit "am schnellsten" erreicht, also nicht abhängig von der Schrittgröße, sondern von der CPU Zeit.
Genau dafür wollte ich die Umrechnung in RA und D haben. Für beobachtende Amateurastronomen ist z.B. eine Genauigkeit von einer Bogensekunde mehr als ausreichend, um die die berechneten Objekte mit Sternkarte und Optik abends am Himmel zu finden. Man hätte dann eine nette und hilfreiche Anwendung :) . Bei den helleren Planeten wäre auch eine Genauigkeit von einer Bogenminute ausreichend. Bei Neptun, Pluto und den Asteroiden sollte die Position aber genauer berechnet werden, um das Auffinden und Identifizieren zu erleichtern.
Gruß

Bernhard

Bernhard
28.05.2011, 12:16
Ich denke viel fehlt da nichtmehr, um zu RA und Dec. von der Erde aus, zu kommen.
Hallo Stefan,

ich habe die Umrechnung jetzt zusammengestellt, allerdings ohne die Abhängigkeit vom Beobachtungsort und damit probehalber die Position der Venus zum 31.03.2011 00:00 mit CNewton::Runge-Kutta (Schrittweite = 1000s) berechnet. Die größere Abweichung ergab sich bei der Deklination zu weniger als 14''. Das ist also gar nicht mal so schlecht :) .

Als nächstes müsste dann der Beobachtungsort inklusive Zeitzone berücksichtigt werden. Insbesondere bei der Position des Mondes sollte der Beobachtungsort berücksichtigt werden.
MfG

Bernhard
28.05.2011, 14:29
Also: Mit Hilfe der Ausgabe von RA und D sind die Ergebnisse leichter und sicherer zu interpretieren. Ich habe deswegen jetzt alle 18 Differenzen mal rechnen lassen: Drei Klassen (CNewton, CNewtonSP und CEIH) mit jeweils drei verschiedenen Methoden (Runge-Kutta mit 1000s, BS mit 8640 s und BS mit 86400s). Berechnung der Position der Venus am 31.03.2012 00:00 in RA und D.

Dabei ergeben sich im Vergleich zu Horizons die folgenden Differenzen in Bogensekunden


<table border=1 cellspacing="0" cellpadding="0">
<tr><td>&nbsp;</td><td colspan="2">RK1000</td><td colspan="2">BS8640</td><td colspan="2">BS86400</td></tr>
<tr><td>&nbsp;</td><td>&Delta; RA/''</td><td>&Delta; D/''</td><td>&Delta; RA/''</td><td>&Delta; D/''</td><td>&Delta; RA/''</td><td>&Delta; D/''</td></tr>
<tr><td>CNewton</td><td>0,16</td><td>21,06</td><td>0,16</td><td>21,06</td><td>0,16</td><td>21,07</td></tr>
<tr><td>CNewtonSP</td><td>41,55</td><td>43,99</td></td><td>41,55</td><td>43,99</td></td><td>41,59</td><td>44,12</td></td></tr>
<tr><td>CEIH</td><td>0,16</td><td>21,04</td><td>0,16</td><td>21,04</td><td>0,16</td><td>21,04</td></tr></table>

Im Rahmen der Genauigkeit von 0,01'' macht es also bei CNewton und CEIH keinen Unterschied, ob mit RK1000, BS8640 oder BS86400 (!) gerechnet wird. Die Vergrößerung des absoluten Fehlers im Vergleich zum vorigen Beitrag von 14'' auf 21'' kommt von der Berücksichtigung der Lichtlaufzeit. Warum der Fehler dadurch größer wird, weiß ich momentan aber auch nicht. Ich habe halt über die Geschwindigkeit der Venus die retardierte Position der Venus linear zurück interpoliert.
MfG

Bernhard
29.05.2011, 19:46
Hallo zusammen,

die verbleibenden Fehler im Bereich mehrerer Bogensekunden schiebe ich mittlerweile doch auf die relativistische Periheldrehung. Zum einen ist es nämlich fraglich, ob der EIH-Limes diese ausreichend genau beschreibt (die numerische Auswertung läßt da eher ein Nein vermuten) und zum anderen braucht man die korrekten Bahnelemente, um die Korrekturen der ART korrekt anzuwenden (ich hatte das weiter oben nur Pi mal Daumen abgeschätzt).

Um nun weiter zu kommen kann man noch die PASCAL-Quelltexte von Cartes Du Ciel lesen. Dort findet man viele Tabellen, die vom JPL stammen. Cartes Du Ciel interpoliert dann zwischen den Werten der Tabelle und kann damit die Werte von Horizons praktisch auf die Bogensekunde genau und besser reproduzieren (!).

Man findet im Netz auch die zugehörige Dokumentation. Ein sehr guter Einstieg dazu ist eine google-Suche nach "J. Chapront, 1995, Astron. Atrophys. Sup. Ser., 109, 181."
Dort findet man einen großen Katalog, als pdf, mit Literatur, wie z.B. http://adsabs.harvard.edu/abs/1994A%26A...282..663S

Die Umrechnung der euklidischen Bahnkoordinaten in RA und D kann man übrigens mit Hilfe einer Translation in geozentrische (bezogen auf den Schwerpunkt der Erde) Koordinaten realisieren. Anschließend dreht man um die x-Achse, die ja in Richtung Frühlingspunkt zeigt. Der Drehwinkel entspricht dabei genau der Schiefe der Ekliptik für J2000.0 und beträgt 23° 26′ 21,45″ (s. http://de.wikipedia.org/wiki/Schiefe_der_Ekliptik#Schwankung_der_Erdachse_und_m ittlere_Ekliptikschiefe). Anschließend transformiert man noch die transformierten euklidischen Koordinaten in Kugelkoordinaten.
MfG

Runzelrübe
29.05.2011, 19:52
Mir fällt es gerade wie Schuppen von den Augen. Vor ein paar Tagen hat Stefan in dem abgesplitteten Thema über eine mögliche Methode, Planeten in Systemen mit bekannten Transitplaneten zu entdecken, ein Paper verlinkt (übrigens ein großes Dankeschön an dieser Stelle!). Ich erinnere mich gerade daran, dort etwas über Unterschiede im Transitzeitpunkt im Sekunden- bis Minutenbereich gelesen zu haben bei Störmassen um Erdgröße.
Habt ihr das System mal inklusive der Jupiter- und Saturnmode laufen lassen? Da kommen wir (über den Daumen gepeilt) doch sicherlich gut und gern auf eine Erdmasse. Horizons scheint in dieser Hinsicht die Gravitationsmassen ja am jeweiligen Körper getrennt zu haben.

Bernhard
29.05.2011, 20:43
Habt ihr das System mal inklusive der Jupiter- und Saturnmode laufen lassen? Da kommen wir (über den Daumen gepeilt) doch sicherlich gut und gern auf eine Erdmasse. Horizons scheint in dieser Hinsicht die Gravitationsmassen ja am jeweiligen Körper getrennt zu haben.
Hi Rübe,

die Planetenmonde sind seit einigen Versionen enthalten. Der Tipp kam weiter oben von UMa und hat auch einiges an Genauigkeit gebracht. Musst einfach mal nach oben blättern. UMa hatte dort ein pdf (xx405.pdf oder so ähnlich) verlinkt, das passende Gm_i-Werte enthält.
MfG

Bernhard
01.06.2011, 17:21
Binary_10 (http://www.coding-facility.com/sr/Binary_10.rar)
Hallo SRMeister,

so finde ich das Gesamtprogramm sehr nett. Eine schöne Visualisierung des Sonnensystems.
Mfg

SRMeister
13.07.2011, 00:17
Hey Ho

Mal ein kleines Update was bei mir die letzten Wochen mit dem Programm passiert ist.
Das gesteckte Ziel war einen neuen Integrator zu entwickeln, bei dem das Hauptaugenmerk auf Geschwindigkeit liegt.
Dabei wollte ich aber (aufgrund der weiter vorne besprochenen Einschränkungen) nicht auf die GPU gehen sondern nur die CPU voll ausreizen. Also was ist damit möglich? Wie bekomme ich das Programm noch schneller? Vom C++ Code her haben wir schon so ziemlich alle möglichen Tricks ausgereizt, das Limit hier war erreicht.
Es blieben also noch folgende Varianten: Den SSE Befehlssatz ausnutzen, den jede CPU seit mind. 10 Jahren sowieso besitzt - der aber kaum Beachtung findet. Damit kann man in einem Durchgang nicht eine, sondern 4 Fließkomma-Berechnungen durchführen, außerdem gibt es bestimme Befehle die SQRT und DIV viel schneller ausführen als die normalen (sind dafür nur etwas ungenauer). Es galt einen Algorithmus zu finden, bei dem man quasi 4 Planeten gleichzeitig berechnet - dazu musste sogut wie die ganze Datenstruktur von Linked Lists auf standart Arrays umgestellt werden (aufgrund von Einschränkungen von SSE - das mag es nur 4 floats einzulesen die nebeneinander im RAM liegen - und auch noch an 16 Byte Grenzen ausgerichtet sind - alles hat seinen Preis)
Kurz und Knapp: ich hab da viel nachgedacht, einen mMn super tollen Algorithmus gefunden, dann umgesetzt - nämlich in Assembler - war aber enttäuscht - da meine Assemblerkenntnisse nicht ausreichend sind um dem C++ Compiler Konkurrenz zu machen - das Programm lief gerade mal 30% schneller - viel weniger als ich mir erhofft hatte.
Dann habe ich das Ganze nochmal programmiert diesmal mithilfe von "SSE Intrinisic" das sind C++ Funktionen die einem quasi die Assembler arbeit abnehmen - meinen Recherchen und anfänglichen Überlegungen zufolge aber schlechter als wenn man gleich in ASM schreibt - aber da täuschte ich mich. mit den Intrinsics lief das Programm gleich 60% schneller (benötigte nurnoch 40% der Zeit) als der gleiche C++ Algorithmus.
Ich habe mir mithilfe des Debuggers den generierten Assemblercode angeschaut und muss sagen - der Compiler macht zwar viel misst, aber scheinbar immernoch weniger als ich. Leider versteh ich da aber *noch* nicht alles. Es bleibt also weiterhin Platz für Optimierung auf diesem Weg.
Der nächste Weg ist aber schon geplant - wozu hat man eine Quad Core CPU ? :) Das Problem dabei ist folgendes: das Programm läuft dann in mehreren Threads die aber einen gemeinsamen Speicher teilen - die Planetenpositionen und -Geschwindigkeiten. Nach jedem Zeitschritt müssen diese Daten also auf einen gemeinsamen Stand gebracht werden. Ein einzelner Zeitschritt ist aber so schnell in der Ausführung: ich schätze etwa 200- 400 CPU Zyklen - ein Schreib- plus -Lesezugriff auf den RAM dauert bei mir aber schon etwa 200 Zyklen - das ist also der Flaschenhals. Keine Ahnung also ob das überhaupt möglich ist, das Programm auf diesem Weg zu beschleunigen. Würde ein Zeitschritt 10.000 Zyklen dauern - würden die 200 Zyklen Lesen + Speichern nicht ins Gewicht fallen und ich hätte eine 4x Beschleunigung.
Keine Ahnung ob man diese Problem lösen kann.
Die Tage werde ich den Quellcode veröffentlichen - ich muss ihn erst noch etwas aufbereiten.

Betrachtet diesen Eintrag mehr als Update Nachricht - mit Wissenschaft hat das leider (noch) nichts zu tun :(

Stefan

Bernhard
13.07.2011, 09:22
Betrachtet diesen Eintrag mehr als Update Nachricht - mit Wissenschaft hat das leider (noch) nichts zu tun :(
Hi Stefan,

sehr interessanter Beitrag und eine gute Motivation die genannten Bibliotheken eventuell auch selber mal anzusehen :) .

Zusätzlich würden mich nochmal die Gründe für die Geschwindigkeitsoptimierung interessieren. Wir hatten ja gesehen, dass der EIH-Algo die relativistischen Effekte schlechter als erwartet berücksichtigt mit der Konsequenz, dass man bereits nach einer Simulationszeit von einem Jahr beobachtbare Abweichungen bekommt.

Möchtest Du als nächsten Schritt die Anzahl der Körper erhöhen? Eventuell noch mehr Asteroiden aufnehmen? Ganz allgemein: Wohin soll die "Reise" gehen?
Gruß

SRMeister
13.07.2011, 10:06
ja also die Relativität hab ich aus diesem Integrator komplett rausgenommen. Er arbeitet im Moment auch mit RK4 ist aber denke ich leicht auf BS umzustellen. (Dazu würd ich dich überreden wollen :) )
Naja jedenfalls möchte ich wie du richtig annimmst, die Anzahl der Körper extrem erhöhen und die Simulation für große Zeiträume gangbar machen. Ich möchte gern einen Simulator schaffen in den man zb. 1000 oder 10000 Planetenembryos steckt und vielleicht mehrere hunderttausend Jahre durchlaufen lässt.
Deshalb hatte ich auch zum Thema RA/ DEC Ephemeriden kein Interesse gezeigt, wobei das Programm mit den letzten Integratoren dafür prinzipiell geeignet erscheint, von der Genauigkeit und Geschwindigkeit her.
Wahrscheinlich muss man um mein Ziel zu erreichen sowieso auf BS gehen, extreme Zeitschritte nehmen und adaptive Zeitschritte umsetzen damit die Kollision noch erkannt werden kann. Ein N-log-N Algorithmus ist denke ich ab 100 Körper oder so, auch schneller.

Edit: Für die SSE Angelegenheit, da hat mir dieser Artikel (http://www.codeproject.com/KB/recipes/sseintro.aspx)sehr geholen. Es gibt auch einen guten (http://www.3dbuzz.com/vbforum/showthread.php?104753-HowTo-Inline-Assembly-amp-SSE-Vector-normalization-done-fast!) für Inline-Assembler-SSE aber aus meinen gemachten Erfahrungen heraus ist das Thema *nicht empfehlenswert*

Bernhard
13.07.2011, 10:51
Ich möchte gern einen Simulator schaffen in den man zb. 1000 oder 10000 Planetenembryos steckt und vielleicht mehrere hunderttausend Jahre durchlaufen lässt.
Benutze dazu doch vorerst die letzte Programmversion ohne Assembler. Dort ist BS ja bereits implementiert und man könnte damit austesten, wie lange das Programm mit N=1000 und t=1000 ... hunderttausend Jahre rechnet. Auch die Programmierung der neuen Startbedingungen ist etwas Arbeit.

Man könnte damit vielleicht schon sehen, ob da etwas zusammenklumpt :confused: .

Erst wenn das ausgetestet ist, würde ich weitere Geschwindigkeitsoptimierungen einbauen.
Gruß

SRMeister
13.07.2011, 11:38
Dabei gibt es 2 Probleme: 1. 1000 Körper ist wahrscheinlich für den (N²) Algorithmus zu viel.
2. Mit großen Zeitschritten und wenig Objekten würde das zumindest jetzt funktionieren, das Problem sind die Kollisionen: dazu braucht man Zeitschritte in denen sich der Körper nur einen Teil seines Radius weiterbewegt damit jede Kollision erkannt wird. Lösung: kleine Zeitschritte nehmen (dann kommt man wieder nicht weit) oder adaptive Zeitschritte: Bei Annäherung werden die Zeitschritte immer kleiner; daraus folgt auch für jedes bisherige Verfahren eine höhere Genauigkeit.

Bernhard
13.07.2011, 12:34
Wie wäre es denn mit diesem Ansatz hier: http://relativ-kritisch.net/forum/viewtopic.php?p=39831#39831. Dort betrachtet man einfach ein Gas unter dem Einfluß der Gravitation. Man könnte dort auch noch eine Zustandsgleichung einbauen, indem man die Dichte mit anderen Größen, wie Druck usw. koppelt.

Für die zeitliche Entwicklung dieses Systems kann man die Lösung in Kugelflächenfunktionen für einen begrenzten Satz an Funktionen entwickeln und dann die Dynamik der Koeffizienten untersuchen. Das ist zwar etwas komplizierter im Ansatz, aber ich denke, dass man damit eigentlich weiter kommen müsste.

Man sollte dazu dann am besten ein neues Thema aufmachen.
Gruß

SRMeister
14.07.2011, 23:22
Hey Bernhard, ich muss zugeben ich hab nur Bahnhof verstanden :(
Gase unter dem Einfluss der Gravitation - soviel hab ich verstanden - dies wird in der Forschung angewendet um die Entstehung von Planetenembryos zu simulieren - quasi die Vorstufe von dem was ich wollte. Da entstehen dann Vortices (glaube das sie es so nannten) und durch Reibung und andere Prozesse verklumpt das Gas/Staub der frühen Scheibe schließlich - ab da ist dieses Modell dann kaum noch zu gebrauchen - höchstens um die Migration der Planeten nach Innen durch Reibung mit der Restgas/staub-Scheibe zu erfassen - in einem hybriden Modell.

Ich bin mir aber gänzlich unsicher ob du das überhaupt meintest.

beste Grüße

Bernhard
15.07.2011, 08:42
Ich bin mir aber gänzlich unsicher ob du das überhaupt meintest.
Hallo Stefan,

so war es schon gemeint. Natürlich interessiert mich die Planetenentstehung aus einer Gas- und Staubwolke auch, und im www gibt es ja auch Simulationen dazu. Ich habe deshalb im AC-Forum einfach mal einen möglichst einfachen Ansatz dazu ohne Materialeigenschaften notiert. Das ergab sich da recht zwanglos aus dem zugehörigen Thema. "Dummerweise" ist das Lösen partieller Differentialgleichungssysteme schon die etwas höhere Kunst und man sollte das momentan wohl einfach den Hauptberuflichen überlassen.

Bei den Kugelflächenfunktionen müsste man den Anfangszustand der "Wolke" auch schon numerisch zerlegen und hätte dabei schon einen recht großen Aufwand. Mit maxima (http://maxima.sourceforge.net/) könnte man das aber mal versuchen.

Bei dem N^2-Algorithmus sehe ich momentan auch kein "Land in Sicht", um Gas- und Staubwolken damit vernünftig zu simulieren . Deshalb auch der Vorschlag eines ganz neuen Ansatzes. Ich meine nur, wir (oder besser ich) stehen ja nicht unter Zeitdruck und deswegen kann man doch ruhig auch mal etwas kompliziertere Ansätze ansehen.
Beste Grüße

SRMeister
15.07.2011, 11:00
Bei den Kugelflächenfunktionen müsste man den Anfangszustand der "Wolke" auch schon numerisch zerlegen und hätte dabei schon einen recht großen Aufwand. Mit maxima (http://maxima.sourceforge.net/) könnte man das aber mal versuchen.
Die Kugelflächenfunktionen stellen doch aber nur die Situation an der Oberfläche dar? Wie komme ich gedanklich zu einer komplexen Staubwolke, in der auch noch mehrere große Objekte stören?



Bei dem N^2-Algorithmus sehe ich momentan auch kein "Land in Sicht", um Gas- und Staubwolken damit vernünftig zu simulieren . Deshalb auch der Vorschlag eines ganz neuen Ansatzes. Ich meine nur, wir (oder besser ich) stehen ja nicht unter Zeitdruck und deswegen kann man doch ruhig auch mal etwas kompliziertere Ansätze ansehen.

Das stimmt, Zeitdruck haben wir zum Glück nicht. Momentan sind wir ja auch eher wieder in der Planungsphase. Ich bin ja auch noch unschlüssig, wie es weitergeht.
Angenommen man schreitet auf dem bisherigen Weg voran, so müsste man den N^2 Algo. ersetzen. Dabei können wir uns dann von der Genauigkeit, die wir erreicht hatten, vollkommen verabschieden.
Trotzdem bin ich unsicher, ob sich das Ziel erreichen lässt, deswegen sehr interessiert an Alternativen.

Bernhard
15.07.2011, 18:36
Die Kugelflächenfunktionen stellen doch aber nur die Situation an der Oberfläche dar? Wie komme ich gedanklich zu einer komplexen Staubwolke, in der auch noch mehrere große Objekte stören?
Man kann zusätzlich die Entwicklungskoeffizienten c_lm als Funktionen vom Radius r und der Zeitkoordinate t betrachten. Die Summation über l läßt man dann von l=0 bis l=N laufen, wobei N die numerische Genauigkeit steuert.

Setzt man diese Entwicklung in die vier partiellen Differentialgleichungen ein, bekommt man ein genähertes, endliches gewöhnliches Differentialgleichungssystem für die vier Entwicklungskoeffizienten c_ilm(r,t) mit i=1,2,3,4.
Gruß

Bernhard
16.07.2011, 00:03
Setzt man diese Entwicklung in die vier partiellen Differentialgleichungen ein, bekommt man ein genähertes, endliches gewöhnliches Differentialgleichungssystem für die vier Entwicklungskoeffizienten c_ilm(r,t) mit i=1,2,3,4.
Hallo SRMeister,

ich fürchte ich muss da leider wieder etwas "zurückrudern". Sowas (Zitat) funktioniert nämlich nur bei dem Laplace-Operator. Wir haben hier aber eine Divergenz und den Gradienten als separate Operatoren und die Kugelflächenfunktionen sind keine Eigenfunktionen dieser Operatoren. D.h. man kann mit dieser Zerlegung die partiellen DGL (Differentialgleichung) leider nicht unmittelbar in eine gewöhnliche DGL transformieren.

Meinen Vorschlag muss ich deswegen vorerst wieder zurückziehen. Zur Lösung der zitierten partiellen DGL muss man vermutlich doch anders vorgehen :confused: .

Ich werde da aber weiter suchen, weil ich die Dynamik der zitierten partiellen DGL durchaus ganz gerne untersuchen/kennenlernen würde...
Gruß

Runzelrübe
17.07.2011, 14:32
Für 1000 zufällige Objekte stößt man abseits der Geschwindigkeitseinbußen noch auf ein paar weitere Herausforderungen.

Die größte davon ist die zufällige Verteilung. Diese sollte auf einem Generator basieren, der mehrere unabhängige Variablen je Objekt produzieren kann. Eine einfache Schleife über denselben Randomizer hilft da nämlich nicht mehr. Am Ende entstehen sonst Körpereigenschaften aus Variablen, die einer zusammenhängenden Teilmenge einer Reihe gehorchen.

Ich bin aber zuversichtlich, dass wir hier gemeinsam sowohl eine Parallelisierung als auch eine noch bessere Pipelineausnutzung hinbekommen können. Die richtige Synchronisierung der Datengrundlage ist dabei gar nicht so schwierig, wie man sich das allgemeinhin vorstellt. Das richtige Warten und Zusammenführen ist viel entscheidender.

Selbst das Stoßproblem stellt keine allzu große Schwierigkeit dar, da alle Körper in einen durch die Geschwindigkeit seiner Körper gewichteten Octtree eingeordnet werden können und damit statt n^2 Tests durchschnittlich nur noch 8n Tests gemacht werden müssen. Die Kollisionstests schüttelt man am Ende einfach noch aus dem Ärmel hinterher und sie bremsen dann ungefähr genauso stark wie ein Deutschlandfähnchen am Auto. :)

Ich suche zu den angesprochenen Themen noch etwas raus und stelle die Links dazu hier rein, wenn ich wieder bei mir bin.

Gruß,
Rübe

SRMeister
18.07.2011, 21:22
Für die zufällige Verteilung hab ich mir folgendes überlegt. Ich hatte in der ersten Version des Programms, bevor ich die Daten der echten Planeten besaß, eine Funktion, die Planeten zufällig verteilt erstellt. Die freien Parameter waren die Entfernung zur Sonne und der Winkel zum Koordinatensystem. Die Funktion berechnete dann die tatsächliche Position (trivial) und die Geschwindigkeit (nicht so trivial) die nötig ist, um den Planeten in einer exakt kreisförmigen Umlaufbahn zu halten.
Man müsste diese Funktion nur um ein paar Parameter erweitern: Exzentrizität und Inklination.
Ich denke, so kann man das von dir angesprochene Problem lösen. Der Randomizer von C++ sollte dafür ausreichen.
Folgendes Thema finde ich interessant.


Die richtige Synchronisierung der Datengrundlage ist dabei gar nicht so schwierig, wie man sich das allgemeinhin vorstellt.
Das stimmt zwar. Die neue C++0x Definition bietet ausgereifte Multithreadingfähigkeiten. Auf diesen würde ich gerne das Programm aufbauen, damit es Plattformunabhängig bleibt (zumindest der Kern)


Das richtige Warten und Zusammenführen ist viel entscheidender.
Solange wir hier von der Nutzung von Multicore Prozessoren reden, ja. Die Parallelisierung mithilfe der SSE SIMD Instructions, funktioniert aber leider nur, wenn die Daten im Speicher linear angeordnet sind, was nur bei O(N^2) der Fall ist.


Selbst das Stoßproblem stellt keine allzu große Schwierigkeit dar, da alle Körper in einen durch die Geschwindigkeit seiner Körper gewichteten Octtree eingeordnet werden können und damit statt n^2 Tests durchschnittlich nur noch 8n Tests gemacht werden müssen.
Bitte definiere Stoßproblem. Wie helfen die Geschwindigkeiten dabei, Kollisionen zu finden?
Nur um mein Problem mit den Kollisionen nochmal deutlich zu machen: Wir benötigen große Zeitschritte, um die Simulation für große Zeiträume fit zu machen. Große Zeitschritte bedeuten, dass sich ähnliche Probleme ergeben, wie beim Tunneln von Quanten. Kollisionen können nur zwischen den Zeitschritten erkannt werden, in etwa vergleichbar mit der Dekohärenz in der Quantenmechanik. Aber ich denke du hast das schon verstanden, nur ich deine Antwort nicht :)

Ich befinde mich momentan auf der Suche nach einem praktikablen Algorithmus. In [1] werden (u.a. von P. Hut - Miterfinder des Barnes-Hut-Algos) viele aktuell genutzte Varianten diskutiert. Leider geht er nicht auf Integratoren ein, sondern nur auf die Organisation der Daten in Baumstrukturen. Außerdem klammern die Autoren scheinbar gezielt die Variablen-Zeitschritt-Methoden aus, die wir in diesem Projekt aber benötigen um Kollisionen zu erkennen. Trotzdem ist das Dokument ein must-read. Wo die Reise ungefähr hingehen wird (denke ich) ist in [2] relativ gut verständlich erklärt. Allerdings ist der Algo mittlerweile schon etwas angestaubt, es gibt sicher modernere Varianten dessen.

Ich freue mich auch auf deine Links!

Stefan

[1] Gravitational N-body Simulations, M. Trenti (1), P. Hut (2) (http://search.arxiv.org:8081/paper.jsp?r=0806.3950&qid=13109987266909a_nCnN_2006680385&qs=n-body)
[2] Barnes Hut Algorithmus in Deutsch erklärt (http://beltoforion.de/barnes_hut/barnes_hut_de.html)

SRMeister
18.07.2011, 21:29
ich fürchte ich muss da leider wieder etwas "zurückrudern". Sowas (Zitat) funktioniert nämlich nur bei dem Laplace-Operator. Wir haben hier aber eine Divergenz und den Gradienten als separate Operatoren und die Kugelflächenfunktionen sind keine Eigenfunktionen dieser Operatoren. D.h. man kann mit dieser Zerlegung die partiellen DGL (Differentialgleichung) leider nicht unmittelbar in eine gewöhnliche DGL transformieren.
Ich hab zwar wieder kein Wort verstanden. Aber bitte schau dir in dem Link [1] von eben unbedingt Abschnitt 3.4 an ( am Besten aber das ganze Doc *daumenhoch*)
Ich vermute das kommt dem ziemlich nahe und finde das auch sehr interessant und bedenkenswert.


Ich werde da aber weiter suchen, weil ich die Dynamik der zitierten partiellen DGL durchaus ganz gerne untersuchen/kennenlernen würde...
Gruß
Da sind wir einer Meinung!

grüße an Alle!

Bernhard
18.07.2011, 22:30
Aber bitte schau dir in dem Link [1] von eben unbedingt Abschnitt 3.4 an ( am Besten aber das ganze Doc *daumenhoch*)
Hallo SRMeister,

ich habe den Preprint gerade mal überflogen und würde das "Ding" erst mal im Archix lassen. Vermutlich muss es solche Zusammenstellungen auch geben, um die Gedanken der zugehörigen Autoren wieder in's Lot zu bringen. Der Leser muss dagegen aufpassen, dass das Abendessen dort bleibt, wo es hingehört. Sorry, für die drastischen Worte, aber sie entsprechen momentan meiner ganz persönlichen Wahrnehmung, was nicht heißen muss, dass man solche Preprints grundsätzlich nicht lesen sollte. Ich bitte dennoch darum solche Themen den Assistenten an den Unis zu überlassen. Die werden schließlich dafür bezahlt, sich so "Zeugs" reinzuziehen.

Sehr interessant finde ich dagegen den Barnes-Hut-Algo. Freundlich und übersichtlich dargestellt mit interessanten Ideen und sehr überzeugenden Resultaten. Klares "Thumbs up". Mein Vorschlag: genau so was umsetzen. Zu den Kollisionen der einzelnen Körper könnte ja vielleicht Runzelrübe etwas weiter helfen. Die können wir IMO nicht einfach so unter den Teppich kehren, ohne an Seriosität zu verlieren :) .

So damit habe ich jetzt auch mal meinen Ärger des Tages hier abgelassen. Ich steh zur Zeit leider etwas unter Druck wegen der Übersetzung für TED. Die möchte ich bis zum Wochenende fertig haben und die hat aktuell definitiv Vorrang vor dem N-Körper-Problem.

Dank an Herrn Deiters für das neue Unterforum.
Gruß

Runzelrübe
19.07.2011, 01:49
Wie helfen die Geschwindigkeiten dabei, Kollisionen zu finden?

Es ist bereits spät, aber ich möchte wenigstens noch in ein paar Sätzen erklären, was es damit auf sich hat.

Erster Ansatz) ungewichteter Octtree:

Wir können den gesamten Raum, in dem die Körper umherschwirren sollen, in acht gleiche Teile unterteilen. Jeder Körper wird sich ja in erster Linie innerhalb der Hülle um alle Körper befinden. Wir setzen in jede der drei Raumrichtungen mal ganz stupide zwei parallele Ebenen um die Hülle, so dass sie einen Quader ergeben, der alle Körper einschließt. Nun achteln wir den Quader - eine Halbierung je Dimension. Betrachten wir die Szene, so können wir zwei Zustände je Körper unterscheiden. Ein Körper befindet sich entweder vollständig innerhalb eines Achtels oder aber er schneidet eine der eben gezogenen Grenzen. Schneidet er eine Grenze, bleibt er der schnittfreien BoundingBox zugeordnet, befindet er sich jedoch vollständig innerhalb eines der acht Knoten, wird er diesem zugeordnet. Wir teilen nun jedes der Achtel immer nur dann weiter, wenn sich mehr als ein Objekt in ihm befindet oder bis wir ein anderes, durch uns festgelegtes Kriterium erreichen.

Was haben wir dadurch erreicht? Nunja, in erster Linie schonmal eine sehr schnelle Möglichkeit, die allernächsten Nachbarn zu finden, ohne uns die Nachbarn je Objekt merken zu müssen.

Zweiter Ansatz) Gewichteter Octtree

Nehmen wir den Baum aus dem ersten Ansatz. Wir berechnen ja in jedem Zeitschritt für jeden Körper seinen Gravitationseinfluss auf die umgebenen Objekte. Für jeden Körper ergibt sich dabei eine neue resultierende Geschwindigkeit. Die Körper bewegen sich weiter. Statt nun jeden Körper auf Kollisionen mit anderen zu testen, schauen wir uns die Bewegung innerhalb unseres Bezugssystems an. Bis eben waren die Positionen der Körper Ist-Werte in unserem Baum. Das wird sich nun ändern. Welche Aussage können wir denn bereits vorab treffen? Nun, wenn überhaupt, dann können sich nur diejenigen Sterne treffen, die sich entweder innerhalb derselben Boundingbox befanden oder die ihre bisherige BoundingBox oder die der direkten Eltern schneiden würden. Wir werden nun zwei Dinge tun. Wir legen um jeden Körper eine weitere, lokale BoundingBox, die die Raumanteile in Bewegungsrichtung widerspiegelt, wenn wir einen großen Zeitschritt nehmen. Nur dann, wenn sich die Bereiche zweier Körper überlappen, kann überhaupt eine Kollision innerhalb derselben BoundingBox stattgefunden haben. Ebenfalls nur dann, wenn eine dieser BBoxen ihre umschließende BBox trifft oder schneidet, muss bei den Nachbarn geprüft werden. Und nur, wenn solch ein Risiko identifiziert wird, muss der Zeitschritt runtergesetzt oder anders verfahren werden. Ansonsten hat da nix und niemand etwas getroffen. Und aus Erfahrung landen so ca. 2 bis 3 Objekte in einem Knoten.

Ein sinnvolles Kriterium wäre demnach noch, dass die Ausdehnung der kleinsten BBox so gewählt wird, dass sie nicht kleiner ist als der durchschnittliche Weg des Zeitschritts.

Falls das jetzt doch zuviel Bahnhof gewesen sein sollte, ich schau mir das morgen nochmal an. :)

Hirschi
19.07.2011, 02:37
Huhu :)
Auch wenn ich von den Details der Sache wenig verstehe, möchte ich vielleicht doch einen Tip in Bezug auf Fortführung der Simulation und Statistik geben. Ich gehe auch davon aus, dass hier die Verklumpung von Staub gemeint ist.
Im Groben dreht sich ja vieles (gerade auch Runzelrübes Ansatz) um Wahrscheinlichkeiten.
Dein gewichteter Octtree ist mMn der richtige Ansatz. Warum? Wenn man die Simulation fortführt mit der nächsten Stufe, 2 verschmolzene Körper fangen 1 verbleibenden ein und jede andere weitergerechnete Konstellation, dann muss! man sowieso gewichten. Die Feinheit des Gewichtungsrasters ist nur abhängig von der gewollten Genauigkeit oder verfügbarer Rechenleistung.
Die relative Geschwindigkeit der benachbarten lokalen BoundingBoxes spielt eine große Rolle. Trägheit ist das Stichwort.
Ungewichtet würde es immer nur zu gleichen lokalen Ausgangsbedingungen führen, obwohl sich subglobal und global schon ganz andere Bedingungen (Geschwindigkeit, Wahrscheinlichkeit einer Massendominanz, Wahrscheinlichkeit eines leeren Achtels) gebildet haben müssen.

SRMeister
21.07.2011, 22:32
Hallo Bernhard,


Sorry, für die drastischen Worte, aber sie entsprechen momentan meiner ganz persönlichen Wahrnehmung

OKAY ;)
Ich konnte jedenfalls einiges an wichtigen Infos rausziehen (komme gleich drauf) außerdem, P. Hut, der "Miterfinder" des Barnes-H-As, hat das Preprint geschrieben! Aber das muss ja nichts heissen.


Sehr interessant finde ich dagegen den Barnes-Hut-Algo.
Ja das stimmt. Interessant ist das, informativ und trotzdem verständlich geschrieben.
Allerdings enthält der Artikel keine Hintergrundinformationen, wozu der BHA geeignet ist und wozu nicht.


... die hat aktuell definitiv Vorrang vor dem N-Körper-Problem.
Kein Problem, ich mach auch ganz langsam im Moment (Erstmal Planen)

So jetzt zu dem arxiv Artikel. Dort wird genauestens darauf eingegangen, wozu die verschiedenen Algos brauchbar sind und wozu nicht. Beispielsweise ist der BHA *nicht besonders* gut geeignet für die Simulation der Entstehung von Planeten, sondern eher für Kosmologie, Kollision von Galaxien, Dunkle Materie Verteilung und sowas. Ich unterstelle dem Herrn Hut mal, dass er weis, wovon er redet.
Er sagt noch, dass bei Galaxien usw. Fehler bei Close Encountern entstehen können, die dazu führen, dass überproportional (unrealistisch) viele Sterne aus dem System geschleudert werden, oder sich in Binärsysteme gegenseitig binden. Er gibt dazu einen sehr interessanten Ausweg, der mich irgendwie stark an MOND erinnert hat: Er dämpft die Gravitation auf kurzen Skalen durch Einführung eines zusätzlichen Terms. Was passiert, wenn man das nicht beachtet, kann man direkt sehen, wenn man im 2. Link (hier nochmal (http://beltoforion.de/barnes_hut/barnes_hut_de.html)) sich ganz unten das Video anschaut: viele Sterne werden raus gekickt, ergo fehlt dem BHA dort, dieser Term.
Leider, und das fehlt mir an dem Text, wird nicht auf Integratoren sowie auf Planetenentstehung eingegangen, es wird nur mehrmals gesagt "dieses und jenes ist dafür ungeeignet" aber Alternativen werden nicht genannt.
Deshalb musste ich mich weiter durch den Arxiv Dschungel kämpfen. Viel uninteressantes, wenig brauchbares findet man da. BTW hab ich auch keine Augen dafür ob ein Artikel "gut oder schlecht" geschrieben ist, ich seh' immer nur, ob ich das dort gesagte Wissen gebrauchen kann/verstehe. Mir fehlt da wohl eine Windung dafür ;)

Ein Preprint, auf dem ich derzeit etwas hängengeblieben bin, ist dieser hier. (http://arxiv.org/abs/0810.3491)
Einige wichtige Essenzen daraus:
-der BS-Integrator wird als extrem langsam eingestuft (daraus habe ich erstmal nur geschlossen dass es Methoden geben muss, von denen ich noch keine Kenntniss hatte)
- Die "grobe" Planeten-Flugbahn wird mithilfe von sogenannten "Symplectic Integrators" berechnet. Irgendwie nehmen die analytische Keplerbahnen. Macht wahrscheinlich Sinn, da die vielen kleinen Körper fast "Probekörper" in der fast statistischen Masseverteilung sind. Dann bleiben nurnoch Close Encounters und Kollisionen als Sonderfälle
- EXTREM wichtig sind variable Zeitschritte. Generell braucht ein langsamerer (von der Sonne entfernterer) Körper weniger Genauigkeit, ein sonnennahes Objekt feinere Zeitschritte. Sonst verschwendet man zuviel Rechenleistung. Bei Kollisionen sind die sowieso notwendig.


Hallo Rübe,



Erster Ansatz) ungewichteter Octtree:
...

Soweit alles klar, das entspricht ja weitgehend dem BHA.


Zweiter Ansatz) Gewichteter Octtree
Falls das jetzt doch zuviel Bahnhof gewesen sein sollte, ich schau mir das morgen nochmal an. :)
Nein, dass war diesmal gut verständlich. Das erinnert mich sehr an die Methoden des Raytracing. Von dort könnte man evtl sogar Funktionen abkupfern. Macht nun aber aufjedenfall Sinn und ist absolut berücksichtigt.

Hi Hirschi,

ich geb dir erstmal keine Antwort sondern hoffe dass Runzelrübe dass übernimmt, da ich mir nicht ganz sicher bin ob ich dich richtig verstehe.

Grüße an Alle!

Bernhard
22.07.2011, 08:42
OKAY ;)
Ich konnte jedenfalls einiges an wichtigen Infos rausziehen (komme gleich drauf) außerdem, P. Hut, der "Miterfinder" des Barnes-H-As, hat das Preprint geschrieben!
Hallo SRMeister,

es kam einfach ein bisschen zuviel auf einmal.

Ich werde mich voraussichtlich aus dem Thema hier jetzt sowieso mehr heraushalten als vorher, weil ich die direkte Planetenenstehung nicht zu meinen Kerninteressen rechne. Für das numerische Lösen von partiellen DGL habe ich schon eher ein Ohr und ein Auge offen, aber das ist halt auch ein 'irre' weites Feld, so dass man so etwas nebenher nur ansatzweise betreiben kann.
Gruß

SRMeister
29.07.2011, 00:00
Ich weiß auch grad nicht so richtig wie es weitergehen soll. Ich sehe ehrlich gesagt nicht, dass ich selbst das Projekt angehen kann.(Planetenentstehung)
Ich hab die letzte Zeit erkundet was nötig wäre und festgestellt: es ist zu viel. Scheitern tu' ich an der Mathematik und Physik.

Also von vorne: Was könnte man realistischerweise machen?
- Kosmologie (Kollision von Galaxien, Dunkle Materie Effeke auf Galaxien, Galaxienhaufenbildung .... ) gibts schon
- Sternkarten - gibt es schon einige
- Integration unseres Sonnensystems - gibts schon von der NASA

Nur mal so ein paar Ideen:
- Radial Velocity Simulator
- Transit Simulator
- ...

*???*

Stefan

Bernhard
29.07.2011, 06:07
Also von vorne: Was könnte man realistischerweise machen?

Hallo Stefan,

mein Vorschlag wäre: vor allem nichts überstürzen. Open-Source-Software ist ein Thema in Entwicklung und muss meiner Erfahrung nach von sehr vielen Usern unterstützt werden, um überleben zu können. Linux ist da das beste Beispiel. Das wird vor allem als Gegenpol zu Windows von sehr vielen Freunden und Liebhabern in ihrer Freizeit gepflegt, gewartet, weiterentwickelt und gesponsort.

Man muss also sehen und abwarten, was die Leser haben wollen und was gebraucht wird. Meine öffentlichen Projekte habe ich mittlerweile alle wieder vom Netz genommen, weil das Feedback einfach ungenügend und entäuschend war. Letztlich erspart mir das aber auch eine Menge freudloser Arbeit.

Zusätzlich kann man das Unterforum hier auch dazu nutzen, um generell Tipps und Erfahrungen mit freier und kostenpflichtiger Software zur Astronomie auszutauschen.
Gruß

UMa
29.07.2011, 10:59
Hallo Stefan,

Ich weiß auch grad nicht so richtig wie es weitergehen soll. Ich sehe ehrlich gesagt nicht, dass ich selbst das Projekt angehen kann.(Planetenentstehung)
...
Also von vorne: Was könnte man realistischerweise machen?


die jetzige Version, Newton-Dynamik mit BS-Integrator, könnte die Veränderung der Bahnen der Planeten berechen, die Auftreten würden, falls ein gelber, roter, weißer, brauner Zwerg oder umherfliegender Planet durch das Sonnensystem fliegen würde. Ähnliches wurde schon mal hier diskutiert
http://www.astronews.com/forum/showthread.php?5173-Roman-Skript-Destabilisierung-eines-Planetensystems.
Es wäre sicher interessent in Abhängigkeit vom Masse, Geschwindigkeit und Abstand des durchfliegenden Objektes die Auswirkungen auf die Bahnen der Planeten des Sonnensystem zu untersuchen. Wie verändern sich die Bahnen der Planeten? Werden Planeten ausgeworfen, kollidieren sie, oder gar vom eindringenden Objekt eingefangen?

Man müsste lediglich zum Sonnensystemsimulator eine weitere Masse hinzufügen, welche in vielleicht 1000AE startet, damit die Anfangsauswirkungen nicht so groß sind, aber die Rechenzeit noch im Rahmen bleibt. Dann ist eine geeignete Geschwindigkeit von sagen wir zwischen 5 und 100 km/s zu wählen, wobei die Richtung ungefähr auf die Sonne zeigt.

Grüße UMa

Bernhard
29.07.2011, 14:00
Man müsste lediglich zum Sonnensystemsimulator eine weitere Masse hinzufügen, welche in vielleicht 1000AE startet, damit die Anfangsauswirkungen nicht so groß sind, aber die Rechenzeit noch im Rahmen bleibt. Dann ist eine geeignete Geschwindigkeit von sagen wir zwischen 5 und 100 km/s zu wählen, wobei die Richtung ungefähr auf die Sonne zeigt.
Nett wäre doch ein Youtube-Filmchen mit einem möglichst anschaulichen, bzw. spannenden Szenario. Das wäre dann ein echtes Astronews-Video. Man müsste dazu nur jemanden (ich werde das aber nicht machen; bin momentan viel zu erkältet für Sowas) finden, der den Code für die grafische Darstellung inklusive schwarzem/gestirntem Hintergrund zur Verfügung stellt.
Gruß

Ich
29.07.2011, 22:17
Nett wäre doch ein Youtube-Filmchen mit einem möglichst anschaulichen, bzw. spannenden Szenario.
!uribin ;-)
Visualisierung ist so ne Sache. Ich hab vor Jahren mal was in OpenGL geschrieben, so mit Planeten und dazwischenrumfliegen. Dabei hab ich erst bemerkt, wie leer das All so ist: nichts außer Schwarz und ab und zu ein kleiner Punkt. Schöner wird's, wenn man die Planeten um Faktor zehn oder so aufbläst.
Aber ich bin kein Programmierer, den Code von damals kann sich keiner anschauen. Ich bastle nur, bis das Prinzip bewiesen ist, und dann reicht's mir wieder.

SRMeister
30.07.2011, 01:23
Oh ja hab vor einiger Zeit mich auch mal mit Opengl beschäftigt (NeHe Tutorials gemacht falls die jemand kennt)

Ich denke es gibt mittlerweile unter .NET einfachere Möglichkeiten sowas in 3D umzusetzen, hab da aber keine Erfahrung.

Die Idee ansich find' ich sehr gut. Man könnte auch auf eine realistische grafische Darstellung vorerst verzichten und einfach in einem Diagramm den Verlauf der Bahnparameter der Planeten darstellen.

Interessant wäre dann auch das resultierende System auf langfristige Stabilität zu untersuchen: dafür ist nämlich das bisherige Verfahren ausreichend (Bulirsch-Stoer)

Bernhard
30.07.2011, 06:46
Ich denke es gibt mittlerweile unter .NET einfachere Möglichkeiten sowas in 3D umzusetzen, hab da aber keine Erfahrung.
Im Prinzip würde auch die 2D-Projektion reichen. Man könnte z.B. die bereits bestehende .NET-Oberfläche noch etwas "aufmotzen": Schwarzer Hintergrund und kleine farbige Kügelchen in leicht unterschiedlicher Größe für die Planeten (Sonne:Gelb,Merkur:Grau, Venus:Weiß,Erde:Blau,Mars:Rot,Jupiter,Saturn:Gelb, Asteroiden:Grau,Uranus,Neptun:Grün,Pluto:Grau). An die Sonne und den "Eindringling" würde ich dann noch die Labels Sonne und WZ (Weißer Zwerg) als Text setzen, dass die klar erkennbar bleiben. Eventuell müsste man auch den Beobachter weg- und heranfahren (Zoom). Die Kennzeichnung der Bahnen könnte man eventuell noch weglassen, um alles übersichtlich zu halten.

Daraus könnte man dann einen kleinen Clip machen. Die einzelnen Frames werden mit einem geeigneten Zeitabstand berechnet und dann zusammengesetzt. Man müsste sich dazu noch nach einer freien Filmbearbeitungssoftware im Netz umsehen (ich denke, da werde ich mich mal ein bisschen umsehen. Gerade das Zusammensetzen finde ich interessant.).

OpenGL würde ich weglassen, weil zu viel Arbeit. Auf weitere Zahlen würde ich auch verzichten, weil ich mir gut vorstellen könnte, dass sich das dann niemand mehr ansieht.
Gruß

SRMeister
30.07.2011, 10:25
sehr gut.
wie wärs wenn man die Bahn der Planeten des letzten halben Jahres oder so, als Schweif anzeigt?
Zoom gibts schon (Mausrad).
Denke mal ich finde dieses WE noch Zeit um ein paar Sachen umzusetzen.

Bernhard
30.07.2011, 11:21
wie wärs wenn man die Bahn der Planeten des letzten halben Jahres oder so, als Schweif anzeigt?
Das wäre vermutlich von Vorteil.

Es gibt zum Clip zusätzlich noch eine entscheidende Hürde: Das Programm müsste Einzelbilder exportieren können. Es müsste z.B. pro zehn Tage (o.ä) ein Bild produzieren, um die Bilderserie dann zu einem Clip zusammenzufügen. Welches Format da geeignet wäre, weiß ich momentan leider nicht, aber ich gehe der Sache nach.
Gruß

SRMeister
30.07.2011, 11:50
Ich kann das Programm so machen das es einen Film auf dem Bildschirm abspielt. Es gibt dann mehrere Programme die Videos vom Bildschirm aufnehmen können, ich glaube sogar unter Windows 7 ist eins integriert. Ich denke die Programme können sicher auch nur Teile des Bildschirms aufnehmen.
Evtl. könntest du also in die Richtung schauen, das wäre wohl am einfachsten.

Bernhard
30.07.2011, 12:29
Hallo Stefan,

vermutlich willst du den Bildschirmspeicher per Programm immer wieder neu beschreiben? Meiner Erfahrung nach, flackert so etwas gerne.

Es gäbe als Alternative kostenlose Konverter, wie z.B. JPGVideo. Das macht aus maximal 32700 durchnummerierten JPG-Dateien eine avi-Datei. Man kann dort auch die Anzahl an Frames pro Sekunde einstellen und könnte damit die Empfehlungen von youtube gut einhalten. Gibt es bei .NET nicht eine Export-Funktion, die ein Bild mit einer vorgegebenen Auflösung als jpg auf die Festplatte rausschreibt?
Gruß

UMa
02.08.2011, 16:24
Hallo Stefan, Bernhard, ...

ich hatte nur wenig Zeit, daher die geringe Beteiligung meinerseits.

Zu SSE, damit habe bisher nur die negative Erfahrung gemacht, dass bei mir nur noch einfache statt doppelte Genauigkeit erzielt wurde, obwohl Hard- und Software angeblich SSE mit double Genauigkeit unterstützen sollten.

Zur Geschwindigkeit, bei mir kostet, wenn ich mich richtig entsinne, das Wurzelziehen beim Berechnen der Entfernungen ca. 40% der gesamten Rechenzeit. Daher speichere ich die Entfernungen in einer Matrix und berechne nur etwas weniger als die Hälfte, da die Entferungen symmetrisch sind.

Zur Parallelisierung gibt es mehrere Möglichkeiten:

Erstens kann man natürlich die Berechnung der Kräfte parallelisieren, da dies am längsten dauert. Dies sollte unabhängig vom algorithmus möglich sein. Der Nachteil ist natürlich, dass das nur für eine große Anzahl von Körpern günstig ist, da sonst der Overhead zu groß wird.

Zweitens kann man beim BS die Berechnungen für eine großen Zeitschritt mit den einzelnen Teilzeitschritten unabhängig von einander ausführen, das keine Ergebisse beim anderen benötigt werden.
Falls man den Kern bestimmen kann, auf welchen die jeweilen Teilprobleme gelöst werden kann eine optimale Beschleunigung erzielt werden.
Für 4 Kerne wäre die Aufteilung so:
Kern 0: 3 schritte + 4 Schritte
Kern 1: 2 Schritte + 5 Schritte
Kern 2: 1 Schritt+ 6 Schritte
Kern 3: 7 Schritte
Nach 7 Zeitschritten, d.h. nach 7 Berechnungen der Entfernungsmatrix ist jeder Kern fertig und die Extrapolation (nichtparallelisert, da im Vergleich zur Berechnung der Entfernungen/Kräfte schnell) kann beginnen.

Wenn man nur 3 Kerne hat, oder nutzen möchte damit man während der Rechnung noch was anderen machen kann, kann man die 7 Schritte weglassen und dadurch etwas weniger extrapolieren. Man wird etwas kürzere Zeitschritte nehmen müssen. Bei 2 Kernen legt man das, was bei 4 Kernen 2 tun zusammen.

Bei 4 Kernen hat man nun pro großem Schritt nur 7 Kraftberechnungen pro Kern statt 28. Bei 3 Kernen sind es 7 statt 21 und bei 2 Kernen sind es 14 statt 28, also optimal.

Wenn man den Kern nicht bestimmen kann, auf dem die Teilberechnungen ausgeführt werden, kann man nur hoffen, dass das Betriebssystem beim Verteilen keinen zu großen Blödsinn macht. Evtl. könnte man durch Veränderung der Reihenfolge der Starts der Berechnungen (z.b. 7 1 2 3 6 5 4 oder so) versuchen Einfluss zu nehmen.

Drittens gibt es natürlich noch die "poor man's" Parallelisierung. Wenn man ohnehin das Programm mehrmals mit verschiedenen Startwerten oder Parametern starten will, startet man das Programm mehrfach gleichzeitig mit diesen Werten.

Bei der Simmulation der Planetenentstehung gibt es den Trick, den Radius der Körper (die Sonne aber nicht, oder nicht so stark) zu vergrößern und so durch häufigere Kollisionen die notwendige Zeit zu verkürzen. In wie weit das da aber unrealisitsch wird, weiss ich nicht. Aber zum Testen oder erste, schnelle Ergebnisse könnte man das einsetzen. Vielleicht am Anfang um einen Faktor 100, (bei der Sonne nicht oder maximal Faktor 10). Das dürfte dann fast 10000mal schneller gehen.

Adaptive Zeitschritte sind bei BS sehr einfach und kosten keine Zeit, da mit dem Unterschied zwischen k und k-1 Extrapolationsschritten ein Fehlerschätzer gratis zur Verfügung steht. Ich denke, dass bei Kollisionen oder beinahe Kollisionen adaptive Zeitschritte notwendig sind.
Bei der adaptive Zeitschrittweite könnte man z.B. 3 Fehlerschranken angeben 0<f1<f2<f3.
Ist der Fehler kleiner als f1 wird der Schritt zwar akzeptiert, aber die Zeitschrittweite im nächsten Schritt etwas vergrößert. Zwischen f1 und f2 ist alles ok, schrittweite bleibt. Zwischen f2 und f3 wird ge eben gemachte Schritt zwar akzeptiert, aber die Zeitschrittweite im nächsten Schritt etwas verkleinert. Ist der Fehler größer als f3 muss der eben berechnet Zeitschritt verworfen und mit einer kleineren Schrittweite wiederholt werden (zur not mehrfach, bis der Fehler nicht mehr größer als f3 ist).

Durch den Zusatzterm wird ermöglicht nahe Vorbeigänge ohne Reduktion der Zeitschrittweite zu berechnen, ohne dass man bei der Kraftberechnung fast durch 0 teilt und im nächsten Schritt der Körper unter hohem unphysikalischem Energiegewinn davon schießt. Im Prinzip wird die Punktmasse dabei über ein größeres Gebiet "vermatscht". Das verhindert die grobe Verletzung des Enegieerhltungssatzes in der Simulation, aber der Stoßwinkel bleibt immer noch falsch. Er ist daher nur für statistische Betrachtungen bei Sternhaufen oder Galaxien geeignet, nicht aber wenn es auf individuelle Körper ankommt.

Sogenannte "Symplectic Integrators" werden verwendet, wenn sich die Bahnen gut durch ein Keplerproblem + kleine Störungen beschreiben lassen. Dann wird das Keplerproblem "exakt" gelöst und nur über die kleinen Störungen integiert. Das ist sehr gut, solange die kleinen Störungen wirklich klein sind. Bei näheren Begegnungen muss dann auf einen anderen Integrator z.B. BS umgeschaltet werden, weil es sonst zu ungenau wird.

BS ist dagegen universell einsetzbar und leichter zu implementieren.
Das BS so viel langsamer ist, kann ich mir nicht vorstellen, da beim BS die Zeitschrittweiten (großer Schritt) größer sind, als bei den "Symplectic Integrators".

Grüße UMa

SRMeister
03.08.2011, 01:15
Zu SSE, damit habe bisher nur die negative Erfahrung gemacht, dass bei mir nur noch einfache statt doppelte Genauigkeit erzielt wurde, obwohl Hard- und Software angeblich SSE mit double Genauigkeit unterstützen sollten.
Das kann ich so nicht bestätigen. Es stimmt zwar, das SSE (SSE1) nur single precision macht, aber praktisch jeder moderne PC seit Pentium 4 (also 2000) kann auch SSE2 womit double möglich ist.
Selbst Intel empfiehlt nurnoch SSEx anstatt die Gleitkommaeinheit zu nutzen. Alle modernen Compiler tuen das auch selbstständig ohne etwas aktivieren zu müssen.
Ein Nachteil gibt es aber: Die Fließkommaeinheit (die FPU) rechnet intern immer mit 80 bit genauigkeit (long double) leider tut das die SSE Einheit nicht, die kann intern auch max. nur 64 bit (double).




Zur Geschwindigkeit, bei mir kostet, wenn ich mich richtig entsinne, das Wurzelziehen beim Berechnen der Entfernungen ca. 40% der gesamten Rechenzeit. Daher speichere ich die Entfernungen in einer Matrix und berechne nur etwas weniger als die Hälfte, da die Entferungen symmetrisch sind.
Das machen wir auch. Ich glaube der Tip kam schonmal von Dir und ist von Bernhard schon lange umgesetzt.



Zur Parallelisierung gibt es mehrere Möglichkeiten:

Erstens kann man natürlich die Berechnung der Kräfte parallelisieren, da dies am längsten dauert. Dies sollte unabhängig vom algorithmus möglich sein. Der Nachteil ist natürlich, dass das nur für eine große Anzahl von Körpern günstig ist, da sonst der Overhead zu groß wird.
das geht mit SSE quasi ohne Overhead, da der Compiler sowieso SSE Befehle benutzt, aber eben nur solche die mit einem Operanden rechnen. Diese kann man theoretisch fast komplett durch Befehle mit 4 Operanden ersetzen. Womit theoretisch kein Overhead entsteht. Die praktische Ausbeute bleibt nur deshalb hinterher weil meine Assemblerfähigkeiten mit dem optimierenden Compiler (noch?) nicht mithalten können.

Extremer Overhead entsteht nur bei Multithreading und fällt deshalb wie du sagst für kleine N aus.




Zweitens kann man beim BS die Berechnungen für eine großen Zeitschritt mit den einzelnen Teilzeitschritten unabhängig von einander ausführen, das keine Ergebisse beim anderen benötigt werden.
...

Das ist im Prinzip korrekt, es müssen halt bei der Extrapolation alle Daten an einem Ort sein, was bedeutet die einzelnen Kerne müssen die Daten über den Speicherbus austauschen, was den großen Overhead erzeugt. An diesem Punkt muss man ganz besonders vorsichtig arbeiten, damit die Synchronisation korrekt von Statten geht.



Wenn man den Kern nicht bestimmen kann, auf dem die Teilberechnungen ausgeführt werden, kann man nur hoffen, dass das Betriebssystem beim Verteilen keinen zu großen Blödsinn macht. Evtl. könnte man durch Veränderung der Reihenfolge der Starts der Berechnungen (z.b. 7 1 2 3 6 5 4 oder so) versuchen Einfluss zu nehmen.
Da würde ich mir keine zu großen Sorgen machen, das Betriebssystem verteilt die Threads optimal. Händisch würde man das nicht besser hinbekommen.



Drittens gibt es natürlich noch die "poor man's" Parallelisierung. Wenn man ohnehin das Programm mehrmals mit verschiedenen Startwerten oder Parametern starten will, startet man das Programm mehrfach gleichzeitig mit diesen Werten.

Bei der Simmulation der Planetenentstehung gibt es den Trick, den Radius der Körper (die Sonne aber nicht, oder nicht so stark) zu vergrößern ... Das dürfte dann fast 10000mal schneller gehen.
Das wäre möglich. Man kann auch noch an anderen Stellen Geschwindigkeit gegen Genauigkeit tauschen. Ich hatte beispielsweise festgestellt, dass durch Einsatz von float statt double der Fehler nur minimal zunimmt. Außerdem hat die SSE Einheit zusätzliche Varianten für 1/x (Reziprok) und 1/Sqrt die die Ergebnisse aus einer Lookup-Tabelle interpolieren und somit bspw. statt 60 Zyklen den gleichen Schritt in 5 Zyklen machen können. Natürlich auf Kosten der Genauigkeit.


Ich denke, dass bei Kollisionen oder beinahe Kollisionen adaptive Zeitschritte notwendig sind.
Genau. Allerdings muss man dann alle Körper in diesem kleinen Zeitschritt machen. Das ist das Problem. Bei Vielen Körpern hat man dann quasi fast immer 2 Körper die sich sehr nahe sind und somit läuft das ganze System ständig in der höchsten Auflösung.


Sogenannte "Symplectic Integrators" werden verwendet, wenn sich die Bahnen gut durch ein Keplerproblem + kleine Störungen beschreiben lassen. Dann wird das Keplerproblem "exakt" gelöst und nur über die kleinen Störungen integiert. Das ist sehr gut, solange die kleinen Störungen wirklich klein sind. Bei näheren Begegnungen muss dann auf einen anderen Integrator z.B. BS umgeschaltet werden, weil es sonst zu ungenau wird.
Ist halt mathematisch und programmtechnisch extrem aufwändig. Von Hamiltonian usw. hab ich auch keine Ahnung und sehe auch nicht dass sich das ändern könnte :)


BS ist dagegen universell einsetzbar und leichter zu implementieren.
Das BS so viel langsamer ist, kann ich mir nicht vorstellen, da beim BS die Zeitschrittweiten (großer Schritt) größer sind, als bei den "Symplectic Integrators".

Es gibt mehrere opensource Simulatoren. Einen davon hab ich mir mal angeschaut. Er bietet mehrere Integratoren darunter auch BS. Der Autor sagt aber der ist nur zu Testzwecken zu gebrauchen weil er so langsam wäre.
Leider konnte ich im Source aber nur die Kommentare verstehen ;) da, aus welchen Gründen auch immer, die Profis heute immernoch in Fortran programmieren :(
Ich schätze, würde man einen wie von dir skizzierten, guten BS-Ansatz unter Ausnutzung der ganzen SIMD und Multicore Hardware in C++/Assembler schreiben, wäre der sicher schneller als den ihr "mathematisch ausgereiztes" System was auf Fortran77 kompiliert wird...

Runzelrübe
05.08.2011, 16:20
Kurze Zusammenfassung eines alten Eintrags:



Ich habe mich [...] mit der 3D Darstellung [...] beschäftigt [...]. Könntet ihr eure Ergebnisse [...] in eine CSV Datei gießen [...], so dass ich eure Ergebnisse [...] mal visuell darstellen kann? Danke!

Da ich ein Pathtracing eingebaut habe [...]


Es muss doch nicht alles doppelt und dreifach gebaut werden, wenn es schon vorhanden ist. ;)

Ich habe lange nicht an der Diskussion teilgenommen, wie mir gerade auffällt. Auch bin ich wohl zu schnell vorgeprescht mit dem Angebot, noch Links rauszusuchen. Irgendwie blieb in den letzten Wochen am Ende des Tages so wenig Zeit übrig. :D

Ich bin zwar kein Freund von Open-Source, kann allerdings mehrere Dinge beisteuern, die uns hier bereits weiterhelfen würden und die ich vollständig in .NET entwickelt habe: einen ASP.NET-GIF-Animator (pseudo-live durch Parametrisierung, leider nur Fullframes, wird sehr schnell sehr groß), einen Video-Animator (exportierbar in ein Wunschformat) und einen 3D-Visualizer (ganz ohne OpenGL oder DirectX). Letzteren hatte ich, wie man sehen kann, bereits vor 3 Monaten angeboten. Nur leider war wohl niemand darauf eingegangen, da der Rest des Eintrags ein Folgethema generiert hatte, das kurzfristig den Fokus verschob.

Mir fällt gerade noch ein, dass ich euch einen klitzekleinen Ausblick auf das, was euch erwartet, sogar jetzt schon geben kann. 2005 kam mir in den Sinn, eine meiner ersten 3D-Web-Visualisierungen per Shockwave zu bauen (sollte ursprünglich die Karte eines Weltraum-Spiels werden): http://www.sternkun.de/galaxy/

Die "Sterne" sind anklickbar und es existiert eine Tastaturunterstützung. Ich kann mir durchaus vorstellen, mich wieder aktiver an dem Thema zu beteiligen. Gebt mir noch ein wenig Zeit. Bis dahin. :)

Bernhard
05.08.2011, 21:10
http://www.sternkun.de/galaxy/
Hallo Runzelrübe,

das sieht interessant aus! Zeigt diese 3D-Karte echte Sterne oder sind die eher zufällig verteilt?

Vorschlag für ein neues Thema: Wie wäre es, diese Karte so auszubauen, dass man dort echte Sterne erkunden kann. Im Koordinatenursprung sitzt die Sonne, bzw. Sagittarius A (Zentrum der Milchstraße). Sichtbar sollten dann die bekanntesten Sterne des Himmels sein, wie z.B. Alpha Centauri, Sirius, Prokyon, Deneb, Wega, Atair usw. usw. Ganz besonders schön wären dann noch die galaktischen Objekte der Milchstraße aus dem Messier-Katalog (und die hellsten aus dem NGC-Katalog). Man könnte dann die Milchstraße wirklich in 3D kennenlernen und erkunden.

Wenn die Seite auf eine definierte Datenstruktur zugreift, könnte man diese Struktur (z.B. Text-Datei) schrittweise von verschiedenen Freiwilligen füllen/erstellen lassen. Das Projekt könnte dann mit der Zeit immer weiter wachsen.
Gruß

Bernhard
01.09.2011, 18:15
Dabei gibt es 2 Probleme: 1. 1000 Körper ist wahrscheinlich für den (N²) Algorithmus zu viel.
Hallo Stefan,

der Vollständigkeit halber möchte ich hier noch ein Paper von J.W. Wadsley (http://hpcc.astro.washington.edu/feedback/gasoline.pdf) verlinken. Dort werden meiner Meinung nach sehr interessante Ergebnisse bestehender Simulationen des N-Körper-Problems mit großem N vorgestellt. Nach kurzem Überfliegen entpuppt sich diese Arbeit als Weiterentwicklung des Barnes-Hut-Algorithmus:


Pkdgrav departed significantly from the original N-body tree code design of Barnes and Hut
Gruß

B.

Bernhard
06.08.2012, 14:03
Hallo SRMeister,

trotz intensiven Bemühens schaffe ich es nicht die letzte Version als Konsolenanwendung für Visual Studio C++ 2010 Express herzustellen. Hast Du eventuell den letzten Code noch komplett, so dass man es ohne Fehler und ohne Warnings übersetzen kann?
Gruß

SRMeister
06.08.2012, 15:20
okay, also lass mich folgende Vorschläge machen:
ich installiere mir jetzt selbst Express, werde meine aktuelle Version darauf testen, anpassen und dann hier veröffentlichen. Die aktuelle Version beinhaltet die von dir zuletzt veröffentliche Version soweit ich mich erinnere, der Integrator ist also auf dem aktuellen Stand.

Dann haben wir erstmal wieder den selben Startpunkt.

Ich habe das Projekt übrigens folgendermaßen weitergeführt:

- reduzieren von 3 Programmiersprachen auf 2 ( core in C++, Windows Teil in C#)
- momentan kein Assembler im core
- momentan keine Konsole im core, core wird als .dll kompiliert und der Windows Teil benutzt diese dann.
- verbesserung des Windowsteil

Es ist also folgendermaßen: Momentan ist der core ein eigenständiges Projekt in der gesamten Projektmappe, den man unabhängig und nur in C++ bearbeiten kann, aber die Schnittstelle muss so bleiben wie sie ist, damit das GUI darauf zugreifen kann.

Ich erkläre das später noch genauer.

Es wäre möglich auch wieder eine Konsole einzusetzen, aber es würde eigentlich wenig Vorteile bringen und ich würde vorschlagen, wir bleiben vorerst bei DLL und GUI als getrennte Komponenten. So kann man auch die Kompetenzen besser aufteilen und man kann durch die Schnittstelle getrennt, unabhängiger arbeiten.

Ich gehe später noch genauer darauf ein, aber evtl wäre ein SVN repository ganz gut. Ich muss mal sehen ob man AnkhSVN mit VS-Express benutzen kann.

Später mehr.

Grüße

Bernhard
06.08.2012, 16:08
Hallo SRMeister,

ich habe den Programmkern mit dem Berechnungscode als Konsolenprojekt für Visual Studio C++ 2010 Express inzwischen bei sourceforge.net eingestellt, damit der Code allen Interessierten zur Verfügung steht: http://sourceforge.net/projects/nbodysimulator/.

SRMeister
06.08.2012, 18:36
Hallo Bernhard,
für den geplanten Branch, (den würde ich gern TheiaSim nennen?)

kann ich den Code so nehmen, alle Integratoren sowie deren relativistische Variante, sowie den Ra/Dec teil, beseitigen, nur den normalen BS übernehmen und dass als Grundlage für ein kombiniertes DLL/Konsolenprojekt nehmen, welches entweder mit dem Windows Teil zusammenarbeitet oder eben als Konsole läuft?

Bernhard
06.08.2012, 18:41
kann ich den Code so nehmen, alle Integratoren sowie deren relativistische Variante, sowie den Ra/Dec teil, beseitigen, nur den normalen BS übernehmen und dass als Grundlage für ein kombiniertes DLL/Konsolenprojekt nehmen, welches mit dem Windows Teil zusammenarbeitet?
Ja klar. Der Code ist (zumindest von mir aus) frei verwendbar und darf auch verändert werden.

SRMeister
06.08.2012, 18:48
Ich frage auch deshalb, weil ich damit eben einen neuen Branch machen will, wobei wir uns vorher, also jetzt darüber einigen, was dort als Grundausstattung enthalten ist.
Ich hoffe ja dass du dann dort auch mit dabei bist :) Wobei wir uns gern später entscheiden können wie wir weitermachen.(wegen adaptiv, kollisionserkennung usw.)

Webmaster
01.10.2012, 14:09
Die hier ursprünglich folgenden Beiträge wurden in ein eigenes Thema ausgelagert, siehe: http://www.astronews.com/forum/showthread.php?6341-Planetenbahnen-in-Excel-programmieren

Bernhard
01.10.2012, 16:08
Die hier ursprünglich folgenden Beiträge wurden in ein eigenes Thema ausgelagert
Danke. .

Bernhard
23.12.2012, 11:25
Passend für ruhige Minuten an Feiertagen noch ein Lesetipp passend zum Thema: http://www.clearskyinstitute.com/xephem/index.html

Das Programm XEphem startete vor vielen Jahren als kleines Linux-Tool und hat sich mittlerweile scheinbar zu einem recht weit verbreiteten Programm gemausert. Die Homepage zeigt in beeindruckender Weise, wie weit man mit tabellarisierten Bahnelementen und der Kepler-Gleichung kommt. Für Untersuchungen zur Stabilität des Sonnensystems über Millionen von Jahren ist das Tool selbstverständlich eher schlecht geeignet.

Der Quelltext ist frei verfügbar.