Bernhard
Registriertes Mitglied
"100 Punkte für den Kandidaten".
Zuletzt bearbeitet:
Sonnensystem: Position der Planeten?
- Ersteller SRMeister
- Erstellt am
Bernhard
Registriertes Mitglied
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ß
Zuletzt bearbeitet:
Bernhard
Registriertes Mitglied
Hallo Stefan,
ich verwende aktuell die folgenden Änderungen für die Euler-Integration:
PHP:
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();
}
PHP:
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 );
}
}
Zuletzt bearbeitet:
Bernhard
Registriertes Mitglied
In der Funktion PSystem::AddObject habe ich ebenfalls die Multiplikation mit Tconst entfernt, damit in der Variable Velocity auch wirklich eine Geschwindigkeit gespeichert wird:
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
PHP:
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
Zuletzt bearbeitet:
SRMeister
Registriertes Mitglied
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
Stefan
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
Stefan
Bernhard
Registriertes Mitglied
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
MfG
Bernhard
Registriertes Mitglied
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
Registriertes Mitglied
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.
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.
Zuletzt bearbeitet:
Bernhard
Registriertes Mitglied
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.
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
Registriertes Mitglied
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
Registriertes Mitglied
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.
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.
Zuletzt bearbeitet:
Bernhard
Registriertes Mitglied
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ß
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
Registriertes Mitglied
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
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
Zuletzt bearbeitet:
Bernhard
Registriertes Mitglied
So leuchtet es ein.Daher ist man auf der sichereren Seite, wenn man immer erst die kleinen Zahlen miteinander addiert und in der Größenordnung wächst.
Danke.
SRMeister
Registriertes Mitglied
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
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
Registriertes Mitglied
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.
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:
mach:
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).
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:
mach:
PHP:
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).
Zuletzt bearbeitet:
Runzelrübe
Registriertes Mitglied
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
Und was ist 2012? Auch richtig! Nein, nicht die Maja, menno, denk ans Schaltjahr!
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]
Zuletzt bearbeitet:
Bernhard
Registriertes Mitglied
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
Registriertes Mitglied
Leider habe ich die Implementierung der Shewchuk'schen Prädikate (hier ein Link zur incircle-Implementierung in 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 / CGAL zurückgreifen.
Viel Erfolg! Ich empfehle CGAL! Dort findet ihr auch viele andere Konzepte.