kleiner planetarer Simulator

Kibo

Registriertes Mitglied
Kann es sein, dass bei k1 - 4 sich Die Position der Planeten nicht ändert?
Das Internet gibt ja so einiges her zum RK4 blos das ist ja schon merkwürdig, wenn ich da jetzt einen Denkfehler (seit über einem Monat ) habe, korrigiere mich bitte:

Ich ging/gehe davon aus, dass das berechnen der Prognosekoordinaten auch irgendwas bringen muss, also rechne ich für P1
Gravitation aus = m1a
Daraus dann m1v = altes v +m1a*schrittweite
Daraus dann Prognosekoordinate P1 = altes x + m1v*schrittweite/2

dann nehme ich dieses P1 (koordinate) und bestimme die gravitation, die der Körper erfahren würde wenn er an dieser stelle steht.
Daraus dann wieder m2v = altes v + m2a*schrittweite
Daraus dann Prognosekoordinate P2 = altes x + m2v*schrittweite/2

Dann bestimme ich wieder die Gravitation an dieser neuen Koordinate und so weiter und so fort.
Da auf Grund dieser linearen Abhängigkeit ich die Koordinaten ja aus der Geschwindigkeit und diese aus der Beschleunigung berechnen muss kann ich die ja auch schlecht unabhängig von einander betrachten :confused:
 
Zuletzt bearbeitet:

Bernhard

Registriertes Mitglied
Ich ging/gehe davon aus, dass das berechnen der Prognosekoordinaten auch irgendwas bringen muss, also rechne ich für P1
Gravitation aus = m1a
Daraus dann m1v = altes v +m1a*schrittweite
Daraus dann Prognosekoordinate P1 = altes x + m1a*schrittweite/2
Nein. Das wäre z.T. ja das Euler-Verfahren. Es gilt vielmehr:
m1v = altes v
und Prognosekoordinate
P1 = altes x + m1v*schrittweite/2.
Das ist der geschickt modifizierte Euler.
Und zusätzlich brauchst Du hier eben auch die Prognosegeschwindigkeit:
P1v = altes v + m1a * schrittweite/2.

Schaue Dir doch mal die Einheiten an! m1a ist eine Beschleunigung in m/s². Wenn man eine Beschleunigung mit der Schrittweite (Einheit Sekunden) multipliziert bekommt man m/s. Also eine Geschwindigkeit. Schon deswegen kann Deine Formel

P1 = altes x + m1a*schrittweite/2
nicht stimmen.

dann nehme ich dieses P1 (koordinate) und bestimme die gravitation, die der Körper erfahren würde wenn er an dieser stelle steht.
Prinzipiell richtig. Vergiss aber auch hier nicht die Prognosegeschwindigkeit.

Daraus dann wieder m2v = altes v + m2a*schrittweite
Daraus dann Prognosekoordinate P2 = altes x + m2a*schrittweite/2
Auch hier wieder der gleiche Fehler, den man schon aufgrund der physikalischen Einheiten aufspüren kann. Es gilt vielmehr:
m2v = P1v
m2a = Gravitation, bzw. Beschleunigung an der Stelle des Prognosepunktes P1.
 

Kibo

Registriertes Mitglied
Morgen Bernhard,

Ich habe deine Änderungen vorgenommen. Leider sind die Spiralen immernoch da, hat sich nicht verbessert.:( Hier mal der gesamte RK:

m1a:=Grav(%Planeten%.x,%Planeten%.y,%Planeten%.z,%PlanetenB%)

%Planeten%.P1v := %Planeten%.h + m1a *Schrittweite (ohne /2 das sonst Bahnen elliptisch werden)

%Planeten%.m1v := %Planeten%.h

%Planeten%.P1:=%Planeten%.x+ %Planeten%.m1v*Schrittweite/2

m2a:=Grav(%Planeten%.P1x,%Planeten%.P1y,%Planeten%.P1z,%PlanetenB%)

%Planeten%.P2v := %Planeten%.h + m2a *Schrittweite

%Planeten%.m2v:=%Planeten%.P1v

%Planeten%.P2:=%Planeten%.x+ %Planeten%.m2v *Schrittweite/2

m3a:=Grav(%Planeten%.P2x,%Planeten%.P2y,%Planeten%.P2z,%PlanetenB%)

%Planeten%.P3v := %Planeten%.x + m3a *Schrittweite

%Planeten%.m3v:=%Planeten%.P2v


%Planeten%.P3:=%Planeten%.x + %Planeten%.m3v *Schrittweite

m4a:=Grav(%Planeten%.P3x,%Planeten%.P3y,%Planeten%.P3z,%PlanetenB%)

%Planeten%.P4v := %Planeten%.h + m4a*Schrittweite

%Planeten%.v:=(%Planeten%.P1v + 2 * %Planeten%.P2v + 2*%Planeten%.P3v + %Planeten%.P4v)/6

mfg
 
Zuletzt bearbeitet:

UMa

Registriertes Mitglied
Hallo Kibo,

so kann das nicht funktionieren. Wo wird z.B. m2a berechnet? Es sind viel zu viele Fehler drin.
Ich habe mal versucht, deine Bezeichnung ungefähr zu übernehmen. Ort steht für den Vektor der Positionen mit 3*n Komponenten. v für den Geschwindigkeitsvektor mit 3*n Komponenten, wobei n die Anzahl der Planeten ist. Wenn du die Planeten getrennt speicherst, musst du jedesmal wenn P1ort= ... P1v=... m1a=... usw. (also überall, außer bei dt=Schrittweite) steht eine Schleife über alle Planeten machen. Du kannst die Planeten nicht nacheinander berechnen.

Grüße UMa

Aus ort_alt und v_alt mache ort_neu und v_neu.
Grav() berechnet die Gravitationsbeschleunigung, abhängig von Ort und Masse.
Klassisches RKV:

dt=Schrittweite
P1ort=ort_alt
P1v =v_alt
m1a=Grav(P1ort)
P2ort=ort_alt + dt/2*P1v
P2v =v_alt + dt/2*m1a
m2a=Grav(P2ort)
P3ort=ort_alt + dt/2*P2v
P3v =v_alt + dt/2*m2a
m3a=Grav(P3ort)
P4ort=ort_alt + dt*P3v
P4v =v_alt + dt*m3a
m4a=Grav(P4ort)
ort_neu=ort_alt + dt/6*(P1v+2*P2v+2*P3v+P4v)
v_neu =v_alt + dt/6*(m1a+2*m2a+2*m3a+m4a)
 

Kibo

Registriertes Mitglied
Hallo UMa,

ich hab die Zeile nachträglich eingefügt. Besten dank für deine Beteiligung

Du kannst die Planeten nicht nacheinander berechnen.

Du schlägst also vor nach der Berechnung von p2 gleich die Berechnung von p2 für alle anderen Planeten mit zu machen?

Auf deine weiteren Hinweise kann Ich erst heute Abend eingehen.



mfg
 
Zuletzt bearbeitet:

Kibo

Registriertes Mitglied
Das klingt einleuchtend, würde auch erklären warum der Euler funktioniert und Heun und RK4 nicht. Da setze ich mich heute Abend gleich ran, spannend:)

mgf
 

UMa

Registriertes Mitglied
Hallo Kibo,
Du schlägst also vor nach der Berechnung von p2 gleich die Berechnung von p2 für alle anderen Planeten mit zu machen?
das ist für alle Verfahren, außer dem Eulerverfahren, zwingend erforderlich. Sonst berechnest du in der Zeile m2a=Grav(P2ort) die Gravitation nicht oder falsch. Du brauchst dafür die P2ort Positionen aller Planeten.

Grüße UMa
 

Kibo

Registriertes Mitglied
Hallo Bernhard,

Danke schonmal. Ich schreib gerade große Teile um und probiere da auch andere Sachen aus. Das ändert aber nicht wirklich was an Genauigkeit oder Features sondern ist nur etwas sauberer und aufgeräumter programmiert. Da ich sowieso die ganze Schleife umstellen muss, bietet sich das ja an. Ich melde mich wenn ich fertig bin. Momentan sehe ich keine Baustelle die ich nicht selber schließen kann, sollte sich das ändern melde ich mich schon.;)

mfg
 

Kibo

Registriertes Mitglied
Nabend, kleines Update aber keine guten Neuigkeiten.

Auch das Berechnen von P1v für alle Planeten vor P2v, sowie P2v vor P3v hat keine Verbesserung gebracht. Die Spirale ist nicht ein bisschen schmaler geworden. Sehr schade, ich habe wirklich gedacht das wäre der Fehler. Ich werde dann wohl mal ein paar Schritte komplett per Hand nachrechnen, vielleicht findet sich so eine Erklärung.:( Aber nicht mehr heute

gute Nacht
 

Kibo

Registriertes Mitglied
Hmm,

Wenn ich UMa's Code anwende dann ist v um den Faktor 4 zu schwach.

Gleiche ich das aus mit:

ort_neu=ort_alt + dt/6*(P1v+2*P2v+2*P3v+P4v) * 4

Dann VERGRÖßERT sich der Abstand der Planeten zueinander nach jedem Orbit um einige %...


PS: Heun macht das selbe, nur ist der Effekt wesentlich schwächer
 
Zuletzt bearbeitet:

Bernhard

Registriertes Mitglied
dt=Schrittweite
P1ort=ort_alt
P1v =v_alt
m1a=Grav(P1ort)
P2ort=ort_alt + dt/2*P1v
P2v =v_alt + dt/2*m1a
m2a=Grav(P2ort)
P3ort=ort_alt + dt/2*P2v
P3v =v_alt + dt/2*m2a
m3a=Grav(P3ort)
P4ort=ort_alt + dt*P3v
P4v =v_alt + dt*m3a
m4a=Grav(P4ort)
ort_neu=ort_alt + dt/6*(P1v+2*P2v+2*P3v+P4v)
v_neu =v_alt + dt/6*(m1a+2*m2a+2*m3a+m4a)
Hallo Kibo,

würdest Du bitte für jede Zeile den konkreten Wert in Deinem Programm bei dem ersten Schritt angeben. Also in der Art

dt = ...
P1ort = (..., ..., ...)
P1v = (..., ..., ...)
usw. bis
v_neu = (..., ..., ...)

dann kann ich nachprüfen, ob alles passt.
MfG
 

Kibo

Registriertes Mitglied
dt= 1000
P1ort=(-7.48e10/0.0/0.0)
P1v =(-11.868006/-84200.040767/0.000000)
m1a=(-17.802000/0.003000/-0.000000)
P2ort=(74800000000.000000/-10525000.000000/0.000000)
P2v =(-17.802000/-105249.998000/0.000000)
m2a= (-0.029670/0.000006/-0.000000)
P3ort= (74799998516.500000/-10525000.000000/0.000000)
P3v = (-11.868000/-63149.998000/0.000000)
m3a= (-0.017802/0.000004/-0.000000)
P4ort= (74799997033.000000/-21049999.500000/0.000000)
P4v = (-5.934000/-21049.999000/0.000000)
m4a= (-0.005934/0.000002/-0.000000)
ort_neu= (74799988131.994064/-84200040.766666/0.000000)
v_neu = (-11.868006/-84200.040767/0.000000)

Nett dass du mal rüber schaust, danke
 

Bernhard

Registriertes Mitglied
Der neue Gesamtvektor (Ort + Geschwindigkeit für beide Sonnen) ergibt sich dann zu

y_1 = (-7.48e10, 2.105e4, 0.0, 5.9308e-3, 2.105e4, 0.0, 7.48e10, -2.105e4, 0.0, -5.9308e-3, -2.105e4, 0.0)

bei einer Schrittweite von 10.000 Sekunden
Der Vollständigkeit halber und zur Kontrolle:

Ich habe das jetzt mal mit "meinem" Programm gerechnet und komme damit auf

y_1 = (-7.48e10, 2.106e8, 0.0, 5.94e1, 2.106e4, 0.0, 7.48e10, -2.106e8, 0.0, -5.94e1, -2.106e4, 0.0)

Bei einigen Komponenten habe ich die Multiplikation mit der Schrittweite übersehen.
 

Kibo

Registriertes Mitglied
Guten Morgen,

Sorry Bernhard, eine der arithmetischen Vektorfunktionen hat anscheinend die v_alt überschrieben kurz bevor ich am ende des Schrittes alle Variablen abgefragt habe. Nach Korrektur des Fehler in der Funktion hat sich aber auch nichts an dem Verhaltend er Planeten gebessert. Dafür sollten die Ausgabewerte jetzt hoffentlich mal stimmen:

dt= 10000
P1ort=/-7.48e10/0.0/0.0
P1v =/0.0/-2.105e4/0.0
m1a=/-0.005934/-0.000000/-0.000000
P2ort=/74800000000.000000/-105250000.000000/0.000000
P2v =/-29.670000/-21050.000000/0.000000
m2a= /-0.005934/0.000008/-0.000000
P3ort= /74799851650.000000/-105250000.000000/0.000000
P3v = /-29.670000/-21049.960000/0.000000
m3a= /-0.005934/0.000008/-0.000000
P4ort= /74799703300.000000/-210499600.000000/0.000000
P4v = /-59.340000/-21049.920000/0.000000
m4a= /-0.005934/0.000017/-0.000000
ort_neu= 74798813199.406601/-841999354.332800/0.000000
v_neu = /-118.680059/-84199.935433/0.000000

Zu beachten ist, dass ich diesmal mit einer Schrittweite von 10000 gerechnet habe.

mfg
 
Oben