Gauss Algorithmus mit Pivotisierung

Neue Frage »

donmemento Auf diesen Beitrag antworten »
Gauss Algorithmus mit Pivotisierung
Hi Leute,

ich möchte eine Dreieckszerlegung durchführen, meine Matrix M ist 16x16, relativ dünn besetzt. Hab mir den Algorithmus programmiert aber die Eigenwerte liefern den falschen Wert für die Determinante (berechne ich dann als Produkt der Hauptdiagonalelemente). Die Werte habe ich mit Maple überprüft. Nachdem ich das mittlerweile drei mal neu programmiert habe verwirrt , seid hir meine letzte Hoffnung

Hier ist mal mein Code, sieht jemand den Fehler ?

! Elimination der Elemente unter der Hauptdiagonalen (Dreieckszerlegung)
do k = 1,(n - 1)
! Pivotsuche
zet = 0
DO i = k,n
IF (abs(M(i,k)) > zet) THEN
zet = abs(M(i,k))
s = i
END IF
END DO
p(k) = s

! Vertauschung
IF (k < s) THEN
DO j = 1,n
zet = M(k,j)
M(k,j) = M(s,j)
M(s,j) = zet
END DO
END IF

! Berechnung der lik
DO i = (k + 1),n
M(i,k) = M(i,k)/M(k,k)
END DO

! Spaltenweise Berechnung von aufdatiertem A

DO j = (k + 1),n
DO i = (k + 1),n
M(i,j) = M(i,j) - M(i,k)*M(k,j)
END DO
END DO
END DO

!Berechnen der Determinante

det = 1

DO j = 1,n
det = det*M(j,j)
END DO

Schon mal vielen Dank für eure Hilfe
Formelmonster Auf diesen Beitrag antworten »

Wenn ich das richtig sehe hast Du vergessen Vorzeichen der Determinante zu beachten:

donmemento Auf diesen Beitrag antworten »

Das ist auf jeden Fall korrekt, allerdings ändert das ja nur das Vorzeichen der Determinante. Wenn ich den Eigenwert suche, für den Det=0 ist, macht das allerdings nichts aus...
tigerbine Auf diesen Beitrag antworten »

Kannst du deine MAtrix mal matlab konform einstellen?
donmemento Auf diesen Beitrag antworten »

kurze Info zu Matrix, es handelt sich um ein Eigenwertproblem (Stabknickung), die Werte EI1-EI4 (Biegesteifigkeiten), F1, F2 (Federsteifigkeiten), F4 (Drehfedersteifigkeit), l1-l4 (Längen) sind Eingabewerte. lambda ist der Eigenwert, für den die Det=0 werden soll...

[0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0;

0,-(EI2*cos(lambda*l2)*lambda^3)/(EI3),lambda^3,0,0,(EI2*sin(lambda*l2)*lambda^3)/(EI3),0,0,0,0,0,0,0,0,0,0;

(EI1*sin(lambda*l1)*lambda^2)/(EI2),0,0,0,(EI1*cos(lambda*l1)*lambda^2)/(EI2),-lambda^2,0,0,0,0,0,0,0,0,0,0;

0,(EI2*sin(lambda*l2)*lambda^2)/(EI3),0,0,0,(EI2*cos(lambda*l2)*lambda^2)/(EI3),-lambda^2,0,0,0,0,0,0,0,0,0;

0,0,(EI3*sin(lambda*l3)*lambda^2)/(EI4),0,0,0,(EI3*cos(lambda*l3)*lambda^2)/(EI4),-lambda^2,0,0,0,0,0,0,0,0;

0,0,0,-sin(lambda*l4)*lambda^2,0,0,0,-cos(lambda*l4)*lambda^2,0,0,0,0,0,0,0,0;

(F4*lambda)/(EI1),0,0,0,-lambda^2,0,0,0,0,0,0,0,(F4)/(EI1),0,0,0;

-(EI1*cos(lambda*l1)*lambda^3)/(EI2),lambda^3,0,0,(EI1*sin(lambda*l1)*lambda^3)/(EI2),-(F1)/(EI2),0,0,0,-(F1)/EI2),0,0,0,0,0,0;

0,0,-(EI3*cos(lambda*l3)*lambda^3)/(EI4),lambda^3,0,0,(EI3*sin(lambda*l3)*lambda^3)/(EI4),-(F3)/(EI4),0,0,0,-(F3)/(EI4),0,0,0,0;

0,0,0,sin(lambda*l4),0,0,0,cos(lambda*l4),0,0,0,1,0,0,0,l4;

-cos(lambda*l1)*lambda,lambda,0,0,sin(lambda l1)*lambda,0,0,0,0,0,0,0,-1,1,0,0;

0,-cos(lambda*l2)*lambda,lambda,0,0,sin(lambda*l2)*lambda,0,0,0,0,0,0,0,-1,1,0;

0,0,-cos(lambda*l3)*lambda,lambda,0,0,sin(lambda*l3)*lambda,0,0,0,0,0,0,0,-1,1;

-sin(lambda*l1),0,0,0,-cos(lambda*l1),1,0,0,-1,1,0,0,-l1,0,0,0;

0,-sin(lambda*l2),0,0,0,-cos(lambda*l2),1,0,0,-1,1,0,0,-l2,0,0;

0,0,-sin(lambda*l3),0,0,0,-cos(lambda*l3),1,0,0,-1,1,0,0,-l3,0]
tigerbine Auf diesen Beitrag antworten »

Ohne Konkrete Zahlen nützt mir das in Matlab natürlich nichts... Augenzwinkern Wollte das zum Vergleich ausrechnen.
 
 
Formelmonster Auf diesen Beitrag antworten »

Nur mal so eine prinzipielle Frage: Musst Du diese Herangehensweise benutzen um die Eigenwerte zu bestimmen? Da Du im wesentlichen Nullstellen eines Polynoms zehnten Grades suchst läuft das ganze auf ein iteratives Verfahren raus. Dieses müsstest Du mindestens zehn mal mit verschiedenen – gut geratenen – Startwerten starten um alle Nullstellen zu finden. In jedem einzelnen Iterationsschritt müsstest Du die Determinante für gegebenes Lambda neu berechnen. Das ist ziemlich teuer!

Wenn Du einfach nur einen beliebigen Algorithmus zur Eigenwertberechnung brauchst solltest Du Dich mal mit dem „QR-Verfahren“ auseinandersetzen.
donmemento Auf diesen Beitrag antworten »

@tigerbine: werte sind frei wählbar, ich hatte ei1-ei4 = 2,1 *10^6, die steifigkeiten alle 0, dann ist der EW = pi/(l1+l2+l3+l4), sind die einzelnen Eulerfälle (man kann die drehfedersteifigkeit -> unendl. gehen lassen für den nächsten euler...

@formelmonster: ich hab das so programmiert, dass ich mit lambda=0 anfange und die det berechne, anschließend lambda=lambda+1, wenn sich das Vorzeichen der Det ändert weiß ich, dass ich an der Nullstelle vorbei bin, dann gehe ich einen schritt zurück mache mit der nächst kleineren zehnerpotenz weiter. So kann ich mich mit beliebiger genauigkeit dem kleinsten EW nähern (ist auch der einzige den ich brauche)

werd mir das QR-Verfahren auf jeden Fall mal ansehen smile
donmemento Auf diesen Beitrag antworten »

erstmal danke für eure Mithilfe.

Bin mittlerweile fest überzeugt, dass der Fehler nicht im Algorithmus, sondern in meiner Matrix ist.
tigerbine Auf diesen Beitrag antworten »

Eine LR-Zerlegung hatte ich hier mal programmiert. Daher nach matlab gefragt. Da du aber Parameter hast, konnte ich das nicht testen. EW Programme habe ich aktuell keine. [WS] - Programmsammlung Numerik Für Gegencheck fehlt mir die Zeit. Da sieht es bei Formelmonster vielleicht besser aus.

Gruß und viel Erfolg. Wink
Neue Frage »
Antworten »



Verwandte Themen

Die Beliebtesten »
Die Größten »
Die Neuesten »