Verfasst am: 03.08.2012, 11:41
Titel: Hilfe zum Thema Dreikörperproblem
Hallo, ich bin im Rahmen eines Projekts zum Thema "Dreikörperproblem" auf einige Probleme gestoßen. Ich bin noch ziemlich neu auf diesem Gebiet deshalb wollte ich euch mal fragen, ob ihr mir nicht helfen könnt.
Code:
G = 6.67384e-11; % Gravitationskonstante
% Masse der beiden Sonnen, Masse des in Bewegung liegenden Planeten kann vernachlässigt werden
m1 = 10e12;
m2 = 10e9 ;
vX_i = sym(rand(2,1)); % Erstellung von 3 Vektoren die die Position des Körpers im Raum definieren.
vX_j = sym(rand(2,1))% Vektor des Planeten
vX_k = sym(rand(2,1));
vP = sym(rand(2,1)); % Neuer Vektor für den zu bewegenden Planeten
vV = sym(rand(2,1)); % Geschwindigkeitsvektor
vTime = linspace(1,1000,100000); % Zeitvektor für die Berechnungen der Planetenbahn
vP = vX_j;
D = distance(vP, vX_i)% Liefert die Distanz zur ersten Sonne
vA = force(G,m1,vP,D)% Berrechnet die Gravitationskraft der ersten Sonne
D = distance(vP, vX_k)% Das gleiche mit der zweiten Sonne
vA = vA + force(G,m2,vP,D)% Addition der Kräfte ala Superposition
axis([-100100-100100])
vV = vV + cntT.*vA; % Die Geschwindigkeit wird immer neu berrechnet, wie berrechnet man delta(cntT)?
vP = vP + cntT.*vV; % der neue Aufenthaltsort wird auch immer neu berrechnet plot(vP(1), vP(2), 'mo'), hold on % Plotten des Aufenthaltspunktes von dem Planeten
Mein Problem ist es, dass die Planetenbahn nie richtig berechnet wird, es liegt vielleicht auch daran, dass ich immer Delta(cntT) brauche aber wie berrechne ich das, weil wenn ich anstatt cntT "2" eingebe, kommt auch nichts sinnvolles dabei raus. Ich wäre auch sehr dankbar, wenn ihr euch das mal ansehen könntet. Habe alles allein gemacht, deswegen glaube ich, dass ich den Fehler nicht sehe.
was mir auffällt, ist die Distanzfunktion. Ich denke, das ist nicht korrekt. Ich würde hier einfach D=norm(v1-v2) verwenden. Es handelt sich dabei um eine Euklidische Norm/Standardnorm?
Du kannst natürlich auch die Distanzfunktion korrigieren:
Wenn ich die Simulation starte, schleudert es den "Planeten" sofort aus dem System. Ich weiß nicht ob es vielleicht an der falschen Datenfütterung liegt. Es sollte ca so aussehen Datei unten angehängt
Vielleicht liegt es ja an einem falschen Vorzeichen oder so.
So, ich habe das ganze jetzt erstmal korrekt mathematisch modelliert. Das ist wesentlich komplizierter als gedacht. Man muss in jedem Fall ein Differentialgleichungssystem der Ordnung 2 lösen. Ich werde das wahrscheinlich heute noch programmieren und posten.
% Masse der beiden Sonnen, Masse des in Bewegung liegenden Planeten kann vernachlässigt werden
m1 = 10e10;
m2 = 10e9 ;
%% Ortsvektoren definieren
% Anfangsort des sich bewegenden Planeten
x0 = rand(2,1);
% Orte der beiden ruhenden Sonnen
x_s1 = rand(2,1);
x_s2 = rand(2,1);
% Betrag der Anfangsgeschwindigkeit
v0_abs=10E-4;
% Geschwindigkeitsvektor des Planeten - Startvektor
v0 = v0_abs.*1/norm((x_s1-x0)+(x_s2-x0))*(1.5*(x_s1-x0)+(x_s2-x0));
% Zeitvektor für die Berechnungen der Planetenbahn % vTime = linspace(0,20,100);
%% Berechnung der Planetenbahn
fig=figure;
hold on;
grid on;
% Plotten der zwei Sonnen
fig=plot(x_s1(1), x_s1(2), 'r*');
fig=plot(x_s2(1), x_s2(2), 'r*');
% Plotten des Anfangsortes des Planeten
fig=plot(x0(1), x0(2), 'ro');
Wieso nicht die Werte, welche ich durch das Integrieren bekomme. von (x0')
Das ist eine gute Frage. Das hatte mich auch erst stutzig gemacht, bis ich begriff, dass es so sein muss. Ich setzt mal ein wenig weiter vorne an. Bei meiner Recherche zum Lösen des Problems, bin ich drauf gestoßen, dass ode45 ausschließlich Differentialgleichungen der Ordnung 1 lösen kann. Daher musste ich zu einem Trick greifen und das ganz transformieren. Die ersten beiden Gleichungen sind daher Hilfsgleichungen, wiederum aber die entscheidenden Gleichungen, denn diese sind zweifach differenziert ausgehend von der Beschleunigung gerade der Ort (xy-Koordinaten) des Planeten. Es gilt x''=a. Daher wird in den Anfangswerten zuerst x0 und dann v0 gesetzt. Wenn wir nun die Rückgabewerte des Solvers betrachten, stehen zwangsläufig die Ortskoordinaten in den ersten beiden Komponenten. Die Geschwindigkeiten in den letzten beiden. Falls du dich dafür interessierst, kannst du diese auch auslesen.
benjo22 hat Folgendes geschrieben:
Oder ist das so gemeint, dass sich die werte von x0 verändern sodass v0 sich dann bei dieser Gleichung wieder verändert
Dieses Konstrukt ist ein Hilfsvektor, der dafür sorgt, dass der Planet sich stets in Richtung der beiden Sonnen bewegt. Es ist das gewichtete Mittel der beiden Richtungsvektoren, normiert und dann mit dem gewünschten Betrag multipliziert.
Wie errechnet man den Wert 1.5. Hab jetzt nämlich das Programm umgeschrieben auf ein Mehrkörpersystem und mir ist aufgefallen, dass die Flugbahn sehr stark von diesem Wert abhängt. Würde mich sehr interessieren sodass ich es jetzt abgleichen kann
Achso, den kannst du weglassen. Ich habe die Vektorsumme damit gewichtet, damit der Planet nicht immer so langweilig genau in das Zentrum der beiden Sonnen "fliegt".
Verschoben: 07.08.2012, 08:34 Uhr von denny Von Objektorientierte Programmierung nach Programmierung
Einstellungen und Berechtigungen
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
RSS