Sonnensystem: Position der Planeten?

[Bernhard]

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]

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

 

[Bernhard]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

Jup, hier haben wir die aktuelle Version.