Verfasst am: 07.12.2016, 23:11
Titel: Verschiedene Ergebnisse einer m-file
Hallo zusammen,
ich bin im Rahmen meiner Bachelorarbeit auf ein dringendes Problem gestoßen.
Ich habe eine m-file entwickelt, welche auf meinem PC ein unterschiedliches Ergebnis bringt wie auf anderen Rechnern. Grundlegend befasse ich mich mit Strukturdynamik und arbeite an einer Methode zur Modellreduktion. Dabei werden Eigenwerte approximiert. Auf meinem PC (Matlab Version R2016a, 64bit Windows 7 Professional 2009) läuft das super ab und ich erhalte eine gute Approximation. Führe ich nun dieselbe m-file auf einem anderen PC aus, so erhalte ich eine deutlich schlechtere Approximation, also ein anderes Ergebnis. Die Version kann nicht das Problem sein, da auf einem anderen Rechner mit R2016a ebenfalls die schlechtere Approximation raus kommt.
Die Struktur beinhaltet Starrkörpermoden, welche u.a. durch den Befehl eig(...) berechnet werden. Ich habe festgestellt, dass bereits bei der Ausführung dieses Befehls eig(...) für die 3 Starrkörpermoden im unten genannten Code leicht abweichende Werte raus kommen, welche meiner Meinung nach durch das weitere Handwerken mit diesen zu dem abweichenden Ergebnis führt. Ich habe diese im Anhang abgelegt.
Ich kann mir die unterschiedlichen Ergebnisse durch verschiedene Rechner nicht wirklich erklären. Vielleicht hätte jemand von euch einen Ansatzpunkt?
Vielen Dank im Voraus für die Hilfe!
hier zur Vollständigkeit die m-file:
Code:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% define all the geometric quantities and material properties
%% Berechnung der free interface normal modes - homogenes Eigenwertproblem
n1=20; %Anzahl der berücksichtigten Eigenmoden in dynamischer Lösung von Struktur 1
n2=15; %Anzahl der berücksichtigten Eigenmoden in dynamischer Lösung von Struktur 2
[Theta_1,lambda_1diag] = eigs(full(B1),-full(A1),n1,'sm');
[lambda_1,k1_EW] = sort(diag(lambda_1diag));
%Zuordnen der Eigenvektoren zu den richtigen Eigenwerten
Theta_1 = Theta_1(:,k1_EW);
[Theta_2ganz,lambda_2diag] = eig(full(B2),-full(A2)); %eigs Befehl aufgrund von Singularität nicht möglich
[lambda_2,k2_EW] = sort(diag(lambda_2diag));
%Zuordnen der Eigenvektoren zu den richtigen Eigenwerten
Theta_2ganz = Theta_2ganz(:,k2_EW);
%% Normalisieren der Eigenvektoren bezüglich A % für komplette Eigenvektoren Struktur 1 [ Theta_1orth ] = Orthonormalize( Theta_1, A1, B1, lambda_1 );
% für komplette Eigenvektoren Struktur 2 [ Theta_2orthganz ] = Orthonormalize( Theta_2ganz, A2, B2, lambda_2 );
%% Entwicklung der regulären und generalisierten Starrkörperformen (Substruktur 2) % Struktur 2 besitzt 3 Starrkörpermoden: % der Verschiebungsvektor besitzt eine reguläre(u) % jedoch keine generalisierte(v) rbm, während der Rotationsvektor für die % eine reguläre und eine generalisierte rbm ausführt. % Somit besitzt das System insgesamt r=3 reguläre und generalisierte rbm
r=3; % Gesamtanzahl der Starrkörpermoden
Theta_2orth = Theta_2orthganz(:,r+1:r+n2); % Abspeichern der normalisierten Eigenvektoren zu EW ungleich 0
Theta_2 = Theta_2ganz (:,r+1:r+n2); % Abspeichern der nicht normalisierten Eigenvektoren zu EW ungleich 0
% nachstehend erfolgt die Manipulation des generalisierten Starrkörpermode % damit eine Menge linear unabhängiger Eigenmoden vorliegt
%% Bool'sche Matrix für Substruktur 1 und Substruktur 2 % Die Geschwindigkeiten der Schnittflächenkoordinaten werden aufgrund der allgemeinen Form der Zustandsraumdarstellung % nicht berücksichtigt (siehe Ausführung)
%% statische Teillösung von Struktur 1 durch Aufprägen von Einheitskräften an den Schnittflächen
%Struktur 1 enthält keine Starrkörpermoden, da keine Eigenwerte = 0
%auftreten. Somit kann die Nachgiebigkeitsmatrix über die Inverse von D gebildet werden
%% statische Teillösung von Struktur 2 durch Aufprägen von Einheitskräften an den Schnittflächen
% Wichtig bei der Wahl der fixierten dofs ist, dass diese zum gleichen % Knoten gehören. % Hier 2 rbm für Rotation und 1 rbm für Durchsenkung. % Es wurde der vorletzte Knoten der 2ten Struktur fixiert. Grundsätzlich % egal welcher innere Knoten fixiert wird (nur kein Knoten, der an einem % Dämpfer befestigt ist)
Ich habe mal reingeschaut. Es fehlen mindestens SystemKraker und Orthonormalize um den Code auszuführen. Der Technische Support würde diesen Code auch benötigen.
Es ist wichtig den Code bezgl. der Unterschiede zu debuggen. Grosse numerische Unterschiede entstehen oftmals durch díe ansteigende Auswirkung von vielen kleinen Unterschieden. Das scheint hier schon teilweise getan zu sein, also wäre z.B. nureine MAT Datei mit den Werten die in "eig" verwendet werden sowie der exakte Aufruf des "eig" Befehls nötig.
Zusätzlich wird des Support bestimmt nach der Ausgabe von
Verfasst am: 08.12.2016, 15:41
Titel: Re: Verschiedene Ergebnisse einer m-file
Hallo FreshFlow,
Im Grunde sind alle Probleme "dringend", die hier gepostet werden.
Die Ausgabe von
eig
kann sich je nach Prozessor, Anzahl der Thread und der verwendeten Bibliothek unterscheiden. Vielleicht benutzt ein Rechner die MKL und der andere nicht?
Wenn die Ergebnisse Deines Programms aber von diesen Rundungsfehlern abhängen, ist der Rechenprozess an sich instabil. In wie weit Du das Resultat dann überhaupt als "Lösung" akzeptierst, ist die Frage. Auf jeden Fall ist es wichtig, die Input-Werte zu variieren im Rahmen der Messgenauigkeit des physikalischen Systems, um dann die Änderung der Outputs zu kontrollieren.
Es ist wie bei der Simulation eines Bleistifts, der auf der Spitze steht: Würdest Du hier dem Ergebnis trauen, dass besagt, dass der Stift für immer gerade stehen bleiben wird?
danke für die schnelle Antwort ich versuche das gleich zu testen.
Bei der Eingabe von ver erscheint das:
----------------------------------------------------------------------------------------------------
MATLAB Version: 9.0.0.341360 (R2016a)
MATLAB License Number: *******
Operating System: Microsoft Windows 7 Professional Version 6.1 (Build 7601: Service Pack 1)
Java Version: Java 1.7.0_60-b19 with Oracle Corporation Java HotSpot(TM) 64-Bit Server VM mixed mode
----------------------------------------------------------------------------------------------------
MATLAB Version 9.0 (R2016a)
Simulink Version 8.7 (R2016a)
Control System Toolbox Version 10.0 (R2016a)
Optimization Toolbox Version 7.4 (R2016a)
Stateflow Version 8.7 (R2016a)
Symbolic Math Toolbox Version 7.0 (R2016a)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
version -blas
ans =
Intel(R) Math Kernel Library Version 11.2.3 Product Build 20150413 for Intel(R) 64 architecture applications, CNR branch AVX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>> version -lapack
ans =
Intel(R) Math Kernel Library Version 11.2.3 Product Build 20150413 for Intel(R) 64 architecture applications, CNR branch AVX
Linear Algebra PACKage Version 3.5.0
Nodes_active1 = unique(Elements(1:10,3:4)); % find which nodes really are used in elements % the degrees of freedom fixed on the ground will % be removed later.-->zählt wie % viele Knoten insgesamt (ohne % z.B.Knoten 1 doppelt zu zählen)
Ndof1 = size(Nodes_active1,1)*2; % 2 dof per node
locnod(Nodes_active1,1:2) = reshape([1:Ndof1],2,Ndof1/2)'; %erstellt ene Matrix in der Zeilenweise die dof der einzelnen Knoten stehen
Nele = size(Elements(1:10,: ),1); %Anzahl der Elemente
K1 = sparse(Ndof1,Ndof1);
M1 = sparse(Ndof1,Ndof1);
[Kel,Mel] = bar_2D(node1,node2,matel);
dof=[locnod(Elements(el,3),[12])...
locnod(Elements(el,4),[12])] ; % only the translational
elseiftype==2,
[Kel,Mel] = beam_2D(node1,node2,matel);
dof=[locnod(Elements(el,3),: )...
locnod(Elements(el,4),: )] ; %bestimmt den Ort, wo Kel bzw Mel in K bzw M eingebaut werden müssen
end
Nodes_active2 = (unique(Elements(11:18,3:4))); % find which nodes really are used in element % the degrees of freedom fixed on the ground will % be removed later.-->zählt wie % viele Knoten insgesamt (ohne % z.B.Knoten 1 doppelt zu zählen)
Ndof2 = size(Nodes_active2,1)*2; % 2 dof per node
locnod(Nodes_active2,1:2) = reshape([1:Ndof2],2,Ndof2/2)'; %erstellt ene Matrix in der Zeilenweise die dof der einzelnen Knoten stehen
Nele = size(Elements(11:18,: ),1); %Anzahl der Elemente
K2 = sparse(Ndof2,Ndof2);
M2 = sparse(Ndof2,Ndof2);
[Kel,Mel] = bar_2D(node1,node2,matel);
dof=[locnod(Elements(el,3),[12])...
locnod(Elements(el,4),[12])] ; % only the translational
elseiftype==2,
[Kel,Mel] = beam_2D(node1,node2,matel);
dof=[locnod(Elements(el,3),: )...
locnod(Elements(el,4),: )] ; %bestimmt den Ort, wo Kel bzw Mel in K bzw M eingebaut werden müssen
end
% fix dofs of nodes 1
%a. find the degrees of freedom that are fixed
dof_fix1 = [];
for n=[1],
dof_fix1 = [dof_fix1 locnod(n,: )];
end
Ndof_fix1 = size(dof_fix1,2);
%b. remove these fixed dof from the list of dof and eliminate them in the matrices
dof_rem1 = setdiff([1:Ndof1],dof_fix1);% remaining degrees of freedom
Ndof_rem1=Ndof1-Ndof_fix1;
K1=K1(dof_rem1,dof_rem1); % from here on, K and M are related to the dof_remaining
M1=M1(dof_rem1,dof_rem1);
Vielen Dank für die Hilfe. Ich konnte letztendlich den Code so abändern, dass er robuster wird. Dazu habe ich die numerischen Nullen der Starrkörpermoden auf exakt 0 gesetzt. So habe ich nun auf allen Rechnern dasselbe Ergebnis.
Vielen Dank aber nochmals für die Unterstützung!
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.