Problemstellung

Während bei Stabilitätsproblemen in der Regel nur der kleinste Eigenwert (Knicklast, Beullast, ...) interessiert, sind bei Schwingungsproblemen neben der Grundfrequenz häufig auch einige weitere der kleinsten Eigenfrequenzen von Interesse. Deshalb wird das Iterationsverfahren zur Berechnung des kleinsten Eigenwerts für das allgemeine Matrizeneigenwertproblem mit symmetrischen Matrizen

Allgemeines Matrizeneigenwertproblem mit quadratischen Matrizen A und B

hier auf diese Aufgabenstellung erweitert. Zunächst wird daran erinnert, dass für dieses Problem gilt (vgl. die Seite "Symmetrische Matrizeneigenwertprobleme"):

  • Alle Eigenwerte sind reell (und positiv, wenn A und B auch noch positiv definit sind).
     
  • Es existiert eine Ähnlichkeitstransformation

Ähnlichkeitstransformation für das allgemeine symmetrische Matrizeneigenwertproblem

    mit der Matrix X, deren Spalten die Eigenvektoren des allgemeinen Matrizeneigenwertproblems enthalten, und der Diagonalmatrix Λ, deren Diagonalelemente die Eigenwerte sind (diese Aussage unterscheidet sich von der über das spezielle Matrizeneigenwertrproblem nur dadurch, dass die Eigenvektoren "B-orthonormal" sind, E ist die Einheitsmatrix).

Inverse Simultaniteration

Die Idee der inversen Simultaniteration: Der Iterationsprozess, der für die Berechnung des zum kleinsten Eigenwert gehörenden Eigenvektors führt, wird mit einem Satz von p Vektoren gestartet (p ≤ n, n ist die Zeilen- bzw. Spaltenanzahl der Matrizen A und B). Die Startvektoren werden zu einer Matrix Y0 zusammengefasst. Wenn man nun den inversen Iterationsprozess für einen Vektor so modifiziert, dass er mit dieser Matrix startet (und dann natürlich in jedem Schritt eine Matrix gleichen Formats Yμ erzeugt), würden alle Spalten der Matrix gegen den Vektor des kleinsten Eigenwerts konvergieren.

Abhilfe schafft folgende Strategie: Nach jedem Iterationsschritt werden die Spalten der iterierten Matrix "von links nach rechts orthonormiert" (für das allgemeine Matrizeneigenwertproblem natürlich "B-orthonormiert"). Dann konvergiert die (jeweils nur "B-normierte") erste Spalte wieder gegen den Eigenvektor des kleinsten Eigenwerts, die jeweils von allen Komponenten des ersten Eigenvektors "gereinigte" zweite Spalte nach der "von-Mises-Idee" gegen den Eigenvektor des zweitkleinsten Eigenwerts usw.

Die Iterationsvorschrift lässt sich also folgendermaßen formulieren:

Iterationsvorschrift für die inverse simultane Vektoriteration für das Allgemeine symmetrische Matrizeneigenwertproblem

Das Verfahren konvergiert gegen die Matrix X, deren Spalten die p Eigenvektoren ("B-orthonormiert") zu den p kleinsten Eigenwerten enthalten. Das Problem der Spaltenorthonormierung kann auf unterschiedliche Weise gelöst werden (Orthonormierungsverfahren nach Schmidt, Subspace-Orthonormierung, ...).

Schmidtsches Orthonormierungsverfahren

OrthonormSchmidt03Das Schmidtsche Orthonormierungsverfahren basiert auf der Idee, die Matrix Z in das Produkt

Zerlegung einer Matris in das Produkt aus einer Orthonormalmatrix und einer Rechtsdreiecksmatrix

mit einer Rechtsdreiecksmatrix R so zu zerlegen, dass für die Matrix Y die Beziehung

Bedingung für "B-orthonormierte" Matrix der Eigenvektoren

gilt ("B-orthonormierte" Spalten).

Rechts wird dies mit dem Falkschen Schema für diese Matrizenmultiplikation veranschaulicht, wobei zu beachten ist, dass nur Z vorgegeben ist. Y und R müssen so bestimmt werden, dass Y der oben angegebenen Gleichung ("B-Orthonormierung") genügt.

Es wird angenommen, dass die Spalten 1 bis k−1 der Matrix Y bereits orthonormiert sind. Die k−te Spalte von Z muss sich dann aus dem Produkt der ersten k Spalten von Y mit den k Elementen rik der k-ten Spalte von R ergeben:

Formel für einen Spaltenvektor

Wenn man diese Beziehung von links mit yiTB multipliziert, dann fallen auf der rechten Seite wegen der "B-Orthonormiertheit" der Spalten alle Produkte bis auf das mit yi, das den Wert 1 hat, weg:

Formel für ein Nichtdiagonalelement der Rechtsdreiecksmatrix

gestattet die Berechnung aller rik der k-ten Spalte von R für i < k (das Hauptdiagonalelement fehlt noch, weil yk noch nicht bekannt ist). Mit den nunmehr bekannten rik kann die oben angegebene Beziehung folgendermaßen umgestellt werden:

Formel für einen Spaltenvektor

gestattet die Berechnung einer vorläufigen (zu den übrigen Spalten orthogonalen, aber noch nicht normierten) k-ten Spalte der Matrix Y. Die Forderung der "B-Orthonormiertheit" kann so formuliert werden:

Berücksichtigung der Bedingung der "B-Orthonormiertheit"

Damit kann der noch fehlende Normierungsfaktor nach

Formel für ein Hauptdiagonalelement der Rechtsdreiecksmatrix

berechnet werden und damit die noch fehlende k-te Spalte von Y:

Endgültige Formel für den nächsten Spaltenvektor der orthonormierten Matrix

Weil die p Spalten von Y gegen die Eigenvektoren der kleinsten p Eigenwerte konvergieren, werden die rik für i < k (die Elemente oberhalb der Hauptdiagonalen) mit fortschreitender Konvergenz immer kleiner. Diese können ja als "Säuberungsfaktoren" aufgefasst werden, die jeweils alle Komponenten der vorläufigen Eigenvektoren yi mit i < k aus dem vorläufigen Eigenvektor yk herausfiltern. Das bedeutet, dass R gegen eine Diagonalmatrix konvergiert, deren Hauptdiagonalelemente die reziproken Eigenwerte enthalten (vgl. die hier angegebene Iterationsvorschrift).

Der hier schrittweise entwickelte Algorithmus wird noch einmal zusammengefasst:

Inverse Simultaniteration für das Allgemeine Eigenwertproblem mit symmetrischen Matrizen:

Start mit einer p-spaltigen Matrix Y (p - Anzahl der zu berechnenden kleinsten Eigenwerte), danach bis zur Konvergenz:

Algorithmus für die inverse Simultaniteration für das Allgemeine symmetrische Matrizeneigenwertproblem

Die Berechnung der Elemente der Dreiecksmatrix R und der Spalten der Matrix Y kann nach dem Schmidtschen Orthonormierungsverfahren erfolgen, hier aufgeschrieben für die Elemente der jeweils k-ten Spalte (k = 1 ... p), der Algorithmus startet für die 1. Spalte mit Schritt 3:

Formelsatz für das Schmidtsche Orthonormierungsverfahren

Bei diesem Algorithmus können Bandstrukturen der beiden Matrizen A und B voll ausgenutzt werden.

Dieser Algorithmus ist realisiert in einem C-Programm isiasb_m.c, für das ein Interface (Mex-File) existiert, mit dem es aus Matlab aufgerufen werden kann. Eine Beschreibung, ein Beispiel und die Möglichkeit zum Download über diesen Link.