| Christian_007 |
 Verfasst am: 28. März 2026 00:02 Titel: Simulation Himmelskörper |
|
Hallo zusammen,
kann bitte jemand meine angehängte m-Datei in Matlab ausführen?
Mir geht es dabei nicht um die Physik, sondern um die Ausführungsgeschwindigkeit in Matlab.
Mit Octave habe ich die Situation, dass mit zunehmender Simulationsdauer, diese immer langsamer wird.
Ich simuliere die Bewegung von drei Himmelskörpern im Raum (Sonne, Erde und Mond). Rein rechnerisch ist das überhaupt keine große Sache. Wenn ich mich nicht verzählt habe, sind das gerade mal 36 Gleichungen, wobei einige Vektorgleichungen dabei sind. Nachdem die Vektoren aber alle nur dreidimensional sind, können es maximal 3 x 36 = 108 Gleichungen sein. Ich würde erwarten, dass das in einer hunderstel Sekunde berechnet werden kann.
Die Simulationsdauer ist nur drei Jahre und die Integrationsschrittweite ein Tag. Das sind demzufolge 3 x 365 = 1095 Tage.
Insgesamt müssen also maximal 108 x 1095 = 118.260 Gleichungen berechnet werden.
Eigentlich hatte ich vor, mehrere hundert oder tausend Jahre zu simulieren, aber mein PC geht schon nach drei Jahren Simulationszeit in die Knie! :-(
Ich habe eine 3D-Grafik (Position der drei Himmelskörper im Raum) und aktuell noch sechs Diagramme.
Mit meinem PC brauche ich fürs erste halbe Jahr Simulationszeit "nur" 15 s, für die drei Jahre aber schon etwa 2,5 min. Hier mal handgestoppte Zeiten:
0,5 - 0:15 min
1 - 0:33 min
1,5 - 0:56 min
2 - 1:23 min
2.5 - 1:55 min
3 - 2:29 min
Ich habe das "Gefühl", dass ein Speicher voll läuft, oder vielleicht immer weiter oder neu "allociert" werden muss. Der PC-Arbeitsspeicher ändert sich aber nahezu gar nicht! Ich habe im Skript auch keine Felder oder dgl., die anwachsen. Wenn ich den Workspace bei Simulationsende abspeichere, dann ist dieser gerade mal 5311 Byte groß! Deshalb tippe ich auf den Grafikspeicher oder ein Problem in Octave?
Je weniger Grafiken man ausgibt, desto schneller läuft die Simulation. Kann ein Grafikspeicher auch volllaufen?
Falls das Skript in Matlab auch knapp 2,5 min für die Simulation von 3 Jahren benötigt: Hat jemand eine Idee, wie man das beschleunigen könnte?
In der Arbeit habe ich z. T. 30 Diagramme (Messdaten), die in Echtzeit dargestellt werden, das ist dort überhaupt kein Problem (bei viel älterer Hardware). Mein PC hier zu Hause ist nur wenige Monate alt.
Vielen Dank schon mal für jede Hilfe!
Christian
PS: Wenn man stabile Himmelsbewegungen simulieren will, dann einfach die Zeilen 40 bis 42 auskommentieren. Ich habe die Zeilen eingefügt, dass die Ergebnisse etwas "spannender" aussehen.
Hier meine m-Datei:
| Code: |
clc;
clear;
close all; % Grafiken löschen
%% *************************************************************************
% Vorbereitungen
einmal = 0; % 1. Durchlauf zur Initialisierung
%% ****************************************************
% Achsen definieren
% 3D-Position im Raum [m]
x_min = -150.0e9; x_max = 155.0e9;
y_min = -150.0e9; y_max = 155.0e9;
z_min = -150.0e9; z_max = 155.0e9;
% resultierende Geschwindigkeit eines Himmelskörpers
v_min = 0; v_max = 80000;
% v_min = -1; v_max = 1; % Skalierung, um die Geschwindigkeitsänderung der Sonne darzustellen
%% *************************************************************************
% Startbedingungen: Positionen der Himmelskörper im Raum, Anfangsgeschwindigkeiten, Massen der Himmelskörper, etc.
t_start = 0; % Startzeit [s]
t_ende = 60*60*24*365 * 3; % Endzeit [s] (x Jahre)
t_ende_J = t_ende / 31536000; % Endzeit in Jahre [a]
dt = 60*60*24 / 1; % Integrationsschrittweite h [s] (1 Tag)
P0S = [ 0.0 ; 0.0 ; 0.0]; % Startposition, Sonne [m]
P0E = [149.6e9 ; 0.0 ; 0.0]; % Startposition, Erde [m]
P0M = [384.4e6 ; 0.0 ; 0.0] + P0E; % Startposition, Mond [m]
v0S_Vektor = [0.0 ; 0.0 ; 0.0]; % Anfangsgeschwindigkeit, Sonne [m/s]
v0E_Vektor = [0.0 ; 29.8e3 ; 0.0]; % Anfangsgeschwindigkeit, Erde [m/s] , stabile Umlaufbahn in der Realität
v0E_Vektor = v0E_Vektor ./ 2; % !!!!!!!!!!!!!!!!!!!!!!!!! Anfangsgeschwindigkeit, Erde [m/s] "verstimmen" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% v0E_Vektor = v0E_Vektor + [0.0 ; 0.0 ; 20e3]; % !!!!!!!!!!!!!!!!!!!!!!!!! Anfangsgeschwindigkeit, Erde [m/s] "verstimmen" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
v0E_Vektor = v0E_Vektor + [0.0 ; 0.0 ; 5e3]; % !!!!!!!!!!!!!!!!!!!!!!!!! Anfangsgeschwindigkeit, Erde [m/s] "verstimmen" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
v0M_Vektor = [0.0 ; 1.02e3 ; 0.0] + v0E_Vektor; % Anfangsgeschwindigkeit, Mond [m/s] , stabile Umlaufbahn in der Realität
mS = 1.989e30; % Masse, Sonne [kg]
mE = 5.9722e24; % Masse, Erde [kg]
mM = 7.35e22; % Masse, Mond [kg]
G = 6.6743015e-11; % Gravitationskonstante [m³/(kgs²)]
% Werte setzen für den ersten Berechnungsschritt setzen (auf Anfangswerte setzen)
PS = P0S; % Position, Sonne im Raum [m]
PE = P0E; % Position, Erde im Raum [m]
PM = P0M; % Position, Mond im Raum [m]
vS_Vektor = v0S_Vektor; % Geschwindigkeitsvektor, Sonne [m/s]
vE_Vektor = v0E_Vektor; % Geschwindigkeitsvektor, Erde [m/s]
vM_Vektor = v0M_Vektor; % Geschwindigkeitsvektor, Mond [m/s]
%% *************************************************************************
% Simulation durchführen
for t = t_start : dt : t_ende % Zeitschleife [s]
%% *************************************
% grafische Ausgabe der Startbedingung
t_J = t / 31536000; % aktuelle Zeit in Jahre für die grafische Ausgabe umrechnen
if einmal == 0 % 1. Durchlauf
scr = get(0, "screensize"); % [x y breite höhe][web:1]
figure("position", [scr(1), scr(2)+70, scr(3)-0, scr(4)-170]); % möglichst großes Ausgabefenster
subplot (3, 3, 1); % 1. Feld (oben links)
view (3); % 3D Sicht einschalten
hold on; % mehr als ein Punkt im Diagramm
grid on; % Gitternetzlinien einschalten
h_S = plot3 (P0S (1), P0S (2), P0S (3), 'ro', 'MarkerSize', 30, 'MarkerFaceColor', 'y');
h_E = plot3 (P0E (1), P0E (2), P0E (3), 'bo', 'MarkerSize', 8, 'MarkerFaceColor', 'b');
h_M = plot3 (P0M (1), P0M (2), P0M (3), 'go', 'MarkerSize', 4, 'MarkerFaceColor', 'g');
xlabel ('x-Position [m]'); % Beschriftung x-Achse
ylabel ('y-Position [m]'); % Beschriftung y-Achse
zlabel ('z-Position [m]'); % Beschriftung z-Achse
axis ([x_min x_max y_min y_max z_min z_max]); % Achsen einstellen
title ('Position der Himmelskörper im Raum: Sonne (gelb), Erde (blau), Mond (grün)'); % Überschrift ausgeben
% resultierende Geschwindigkeiten (Skalar) Erde, Sonne, Mond im 2. Diagramm
subplot (3, 3, 2);
hold on;
grid on;
plot (t_J, norm (vS_Vektor), 'c'); % res. Geschwindigkeit, Sonne [m/s]
plot (t_J, norm (vE_Vektor), 'b'); % res. Geschwindigkeit, Erde [m/s]
plot (t_J, norm (vM_Vektor), 'g'); % res. Geschwindigkeit, Mond [m/s]
axis ([t_start t_ende_J v_min v_max]); % Achsen einstellen
xlabel ('Zeit [a]'); % Beschriftung x-Achse (Zeit [s])
ylabel ('Geschwindigkeit [m/s]'); % Beschriftung y-Achse
title ('v-Erde (blau), v-Sonne (türkis), v-Mond (grün)'); % Überschrift ausgeben
% resultierende Kräfte auf Himmelskörper (Skalar) Erde, Sonne, Mond
subplot (3, 3, 3);
hold on;
grid on;
plot (t_J, 0, 'b'); % res. Kraft, Sonne-Erde [N] (Wert noch nicht bekannt)
plot (t_J, 0, 'c'); % res. Kraft, Erde-Mond [N] (Wert noch nicht bekannt)
plot (t_J, 0, 'g'); % res. Kraft, Sonne-Mond [N] (Wert noch nicht bekannt)
axis ([t_start t_ende_J 0 2e22]); % Achsen einstellen
xlabel ('Zeit [a]'); % Beschriftung x-Achse (Zeit [s])
ylabel ('resultierende Kraft [N]'); % Beschriftung y-Achse
title ('Res. Kraft zwischen Himmelskörpern - Sonne-Erde (blau), Erde-Mond (türkis), Sonne-Mond (grün)'); % Überschrift ausgeben
% Position Sonne
subplot (3, 3, 4);
plot (t_J, PS (1), 'r');
plot (t_J, PS (2), 'g');
plot (t_J, PS (3), 'b');
axis ([t_start t_ende_J -4e6 4e6]); % Achsen einstellen
xlabel ('Zeit [a]'); % Beschriftung x-Achse
ylabel ('Position, Sonne [m]'); % Beschriftung y-Achse
title ('Position Sonne - x-Position (rot), y-Position (grün), z-Position (blau)'); % Überschrift ausgeben
grid on;
hold on;
% Position Erde im Raum
subplot (3, 3, 5);
hold on;
grid on;
plot (t_J, PE (1), 'r'); % x-Position der Erde in rot
plot (t_J, PE (2), 'g'); % y-Position der Erde in grün
plot (t_J, PE (3), 'b'); % z-Position der Erde in blau
axis ([t_start t_ende_J x_min x_max]); % Achsen einstellen
xlabel ('Zeit [a]'); % Beschriftung x-Achse
ylabel ('Position, Erde [m]'); % Beschriftung y-Achse
title ('Position Erde [m] - x-Position (rot), y-Position (grün), z-Position (blau)'); % Überschrift ausgeben
% Position Mond
subplot (3, 3, 6);
hold on;
grid on;
plot (t_J, PM (1), 'r'); % x-Position des Mondes in rot
plot (t_J, PM (2), 'g'); % x-Position des Mondes in grün
plot (t_J, PM (3), 'b'); % x-Position des Mondes in blau
axis ([t_start t_ende_J x_min x_max]); % Achsen einstellen
xlabel ('Zeit [a]'); % Beschriftung x-Achse
ylabel ('Position, Mond [m]'); % Beschriftung y-Achse
title ('Position Mond [m] - x-Position (rot), y-Position (grün), z-Position (blau)'); % Überschrift ausgeben
% Entfernungen
subplot (3, 3, 7);
hold on;
grid on;
plot (t_J, 0, 'b'); % x-Position des Mondes in rot
plot (t_J, 0, 'c'); % x-Position des Mondes in grün
plot (t_J, 0, 'g'); % x-Position des Mondes in blau
axis ([t_start t_ende_J 0 x_max]); % Achsen einstellen
xlabel ('Zeit [a]'); % Beschriftung x-Achse
ylabel ('Entfernung [m]'); % Beschriftung y-Achse
title ('Entfernungen zwischen den Himmelskörpern [m] - Sonne-Erde (blau), Erde-Mond (türkis), Sonne-Mond (grün)'); % Überschrift ausgeben
einmal = 1; % einmaligen Durchlauf nicht wiederholen
else
%% *************************************************************************
% Positionsberechnungen
% Gravitationskraft: Sonne-Erde
dP_SE_Vektor = PS - PE; % Abstandsvektor Sonne-Erde
dP_SE_Vektor_quadrat = dP_SE_Vektor.^2; % Quadrat des Abstandes
r_SE = (dP_SE_Vektor_quadrat (1) + dP_SE_Vektor_quadrat (2) + dP_SE_Vektor_quadrat (3))^(1/2); % Abstand (Skalar) berechnen
F_SE = G * mS * mE / r_SE^2; % Grafitationskraft Sonne-Erde (Skalar) [N]
F_SE_Vektor = F_SE .* (dP_SE_Vektor ./ r_SE); % Kraftvektor in x-, y-, und z-Richtung
% Gravitationskraft: Erde-Mond
dP_EM_Vektor = PE - PM; % Abstandsvektor Erde-Mond
dP_EM_Vektor_quadrat = dP_EM_Vektor.^2; % Quadrat des Abstandes
r_EM = (dP_EM_Vektor_quadrat (1) + dP_EM_Vektor_quadrat (2) + dP_EM_Vektor_quadrat (3))^(1/2); % Abstand (Skalar) berechnen
F_EM = G * mE * mM / r_EM^2; % Grafitationskraft Erde-Mond (Skalar) [N]
F_EM_Vektor = F_EM .* (dP_EM_Vektor ./ r_EM); % Kraftvektor in x-, y-, und z-Richtung
% Gravitationskraft: Sonne-Mond
dP_SM_Vektor = PS - PM; % Abstandsvektor Sonne-Mond
dP_SM_Vektor_quadrat = dP_SM_Vektor.^2; % Quadrat des Abstandes
r_SM = (dP_SM_Vektor_quadrat (1) + dP_SM_Vektor_quadrat (2) + dP_SM_Vektor_quadrat (3))^(1/2); % Abstand (Skalar) berechnen
F_SM = G * mS * mM / r_SM^2; % Grafitationskraft Sonne-Mond (Skalar) [N]
F_SM_Vektor = F_SM .* (dP_SM_Vektor ./ r_SM); % Kraftvektor in x-, y-, und z-Richtung
% **** Erde ****
% Beschleunigungsvektor
F_E_Vektor = F_SE_Vektor + F_EM_Vektor; % Resultierende Kraft auf die Erde berechnen
aE_Vektor = F_E_Vektor / mE; % Beschleunigungsvektor berechnen (F = m * a -> a = F / m) [m/s²]
% Geschwindigkeitsänderung (Vektor) berechnen
dvE_Vektor = aE_Vektor * dt; % Geschwindigkeitsänderung: dv = a * dt [m/s]
% Geschwindigkeitsvektor, Erde berechnen
vE_Vektor = vE_Vektor + dvE_Vektor; % Geschwindigkeit: v = v + dv [m/s]
vE = norm (vE_Vektor); % Geschwindigkeit, Erde (Skalar) [m/s]
% Wegänderung (Vektor), Erde berechnen
dxE_Vektor = vE_Vektor .* dt; % Wegänderung: v = s/t -> ds = v * dt [m]
% neue Koordinaten, Erde berechnen
PE = PE + dxE_Vektor; % Positionsänderung: x = x + dx [m]
% **** Sonne ****
% Beschleunigungsvektor
F_S_Vektor = F_SE_Vektor + F_SM_Vektor; % Resultierende Kraft auf die Sonne berechnen
aS_Vektor = F_S_Vektor / mS; % Beschleunigungsvektor berechnen (F = m * a -> a = F / m) [m/s²]
% Geschwindigkeitsänderung (Vektor) berechnen
dvS_Vektor = aS_Vektor * dt; % Geschwindigkeitsänderung: dv = a * dt [m/s]
% Geschwindigkeitsvektor, Erde berechnen
vS_Vektor = vS_Vektor + dvS_Vektor; % Geschwindigkeit: v = v + dv [m/s]
vS = norm (vS_Vektor); % Geschwindigkeit, Sonne (Skalar) [m/s]
% Wegänderung (Vektor), Erde berechnen
dxS_Vektor = vS_Vektor .* dt; % Wegänderung: v = s/t -> ds = v * dt [m]
% neue Koordinaten, Erde berechnen
PS = PS + dxS_Vektor; % Positionsänderung: x = x + dx [m]
% **** Mond
% Beschleunigungsvektor
F_M_Vektor = F_EM_Vektor + F_SM_Vektor; % Resultierende Kraft auf den Mond berechnen
aM_Vektor = F_M_Vektor / mM; % Beschleunigungsvektor berechnen (F = m * a -> a = F / m) [m/s²]
% Geschwindigkeitsänderung (Vektor) berechnen
dvM_Vektor = aM_Vektor * dt; % Geschwindigkeitsänderung: dv = a * dt [m/s]
% Geschwindigkeitsvektor, Erde berechnen
vM_Vektor = vM_Vektor + dvM_Vektor; % Geschwindigkeit: v = v + dv [m/s]
vM = norm (vM_Vektor); % Geschwindigkeit, Mond (Skalar) [m/s]
% Wegänderung (Vektor), Erde berechnen
dxM_Vektor = vM_Vektor .* dt; % Wegänderung: v = s/t -> ds = v * dt [m]
% neue Koordinaten, Erde berechnen
PM = PM + dxM_Vektor; % Positionsänderung: x = x + dx [m]
%% *******************************************************
% grafische Ausgabe der Punkte (Verschiebung, etc.)
subplot (3, 3, 1); % Positionen Himmelskörper im Raum
set (h_S, 'XData', PS (1), 'YData', PS (2), 'ZData', PS (3)); % Position der Sonne aktualisieren
set (h_E, 'XData', PE (1), 'YData', PE (2), 'ZData', PE (3)); % Position der Erde aktualisieren
set (h_M, 'XData', PM (1), 'YData', PM (2), 'ZData', PM (3)); % Position der Sonne aktualisieren
% plot3 (x, y, z, 'b.'); % Spur zeichnen
drawnow; % Grafik aktualisieren
subplot (3, 3, 2); % Geschwindigkeiten Sonne, Erde, Mond
plot (t_J, vS, 'c');
plot (t_J, vE, 'b');
plot (t_J, vM, 'g');
subplot (3, 3, 3); % resultierende Kraft auf Himmelskörper
plot (t_J, F_SE, 'c'); % Sonne-Erde
plot (t_J, F_EM, 'b'); % Erde-Mond
plot (t_J, F_SM, 'g'); % Sonne-Mond
subplot (3, 3, 4); % Position, Sonne
plot (t_J, PS (1), 'r');
plot (t_J, PS (2), 'g');
plot (t_J, PS (3), 'b');
subplot (3, 3, 5); % Position, Erde
plot (t_J, PE (1), 'r');
plot (t_J, PE (2), 'g');
plot (t_J, PE (3), 'b');
subplot (3, 3, 6); % Position, Mond
plot (t_J, PM (1), 'r');
plot (t_J, PM (2), 'g');
plot (t_J, PM (3), 'b');
subplot (3, 3, 7); % Position, Mond
plot (t_J, r_SE, 'b');
plot (t_J, r_EM, 'c');
plot (t_J, r_SM, 'g');
endif
% pause(0.001); % Geschwindigkeit der Animation
end
|
|
|
