MATHEMATISCHES INSTITUT DER UNIVERSITÄT MÜNCHEN Dr. W. Spann, G. Kleen WS 2022/2023

Programmieren II für Studierende der Mathematik

Aufgabe 10

Das zyklische Jacobi-Verfahren dient zur Berechnung der Eigenwerte und Eigenvektoren einer symmetrischen Matrix A= (a_{ij})_{i,j=0,...,n−1} ∈ ℝ^{n×n} .

Dazu wird die Matrix A sukzessive den Transformationen Q^{⊤}AQ unterworfen, worin Q := Q^{ij} (i < j) eine Matrix bezeichnet, die sich von der Einheitsmatrix nur in den folgenden 4 Einträgen unterscheidet:

\begin{matrix}q_{ii} = & q_{jj} = & c, \\ q_{ij} = & −q_{ji} = & s. \end{matrix}

c und s erfüllen c^{2} + s^{2} = 1 und werden für a_{ij}^{alt} ≠ 0 so gewählt, daß a_{ij}^{neu} = 0, c > 0, und s ∈ (−c, c] gilt. Wenn bereits a_{ij}^{alt} = 0 erfüllt ist, setzt man c = 1 und s = 0. Der Ansatz c = st führt dann nach einigen elementaren Umformungen zu

τ = \frac{a_{jj}^{alt} − a_{ii}^{alt}}{2a_{ij}^{alt}} , t = \frac{sign(τ)}{|τ| + \sqrt{1 + τ^{2}}} mit sign(τ) = \{ \begin{matrix}1 & : & τ ≥ 0 \\ −1 & : & τ < 0 \end{matrix} ,

c = \frac{1}{\sqrt{1 + t^{2}}} , s = ct

Dabei werden nur die i-te und j-te Zeile und Spalte von A^{alt} verändert. Bezeichnet man die k-te Zeile bzw. k-te Spalte einer Matrix M mit m^{k} bzw. m_{k}, so ergibt sich mit A′ := A^{alt}Q und A^{neu} = Q^{⊤}A′

\begin{aligned}a′_{i} & = & ca_{i}^{alt} − sa_{j}^{alt} \\ a′_{j} & = & sa_{i}^{alt} + ca_{j}^{alt} \\ (a^{neu})^{i} & = & c(a′)^{i} − s(a′)^{j} \\ (a^{neu})^{j} & = & s(a′)^{i} + c(a′)^{j} \end{aligned}

Ebenso kann der Update des Produkts U aller bisher angewendeten Matrizen Q über U^{neu}=U^{alt}Q erfolgen:

\begin{aligned}u_{i}^{neu} & = & cu_{i}^{alt} − su_{j}^{alt} \\ u_{j}^{neu} & = & su_{i}^{alt} + cu_{j}^{alt} \end{aligned}

Die Summe der Betragsquadrate der Nichtdiagonalelemente, die im Abbruchkriterium benötigt wird, kann folgendermaßen berechnet werden:

N(A^{neu}) = N(A^{alt}) − 2|a_{ij}^{alt}|^{2}

In einem Zyklus werden zeilenweise nacheinander die Indexpaare (i, j) der oberen Hälfte von A durchlaufen (d.h. j = i + 1,…,n − 1 für i = 0,…,n − 1).

Die Berechnung soll abgebrochen werden, wenn N(A^{neu}) := 2 ∑_{i=0}^{n−1} ∑_{j=i+1}^{n−1} |a_{ij}^{neu}|^{2} ≤ ε · N(A) nach

einem Zyklus gilt oder eine vorgebenene Zahl maxzykl von Zyklen überschritten ist.

Die Näherungen für die Eigenwerte stehen in der Diagonale von A^{neu}, die Näherungen für die Eigenvektoren ergeben sich aus den Spalten der Matrix U, die das Produkts aller bisher angewendeten Matrizen Q ist.

Erstellen Sie eine Funktion jacobi, die das zyklische Jacobi-Verfahren implementiert. Ihr Rückgabewert soll über den Erfolg des Verfahrens Auskunft geben. A, ε und maxzykl sollen über Eingabeparameter, die berechneten Eigenwerte als Vektor und die Eigenvektoren als Matrix über Ausgabeparameter übergeben werden.

Die Matrizen sollen einer eigens definierten Klasse Matrix angehören, der der Datentyp valarray<double> zugrunde liegt, und die Vektoroperationen für Matrixzeilen bzw. -spalten zulässt. (Ergänzen Sie z.B. den in der Vorlesung vorgestellten Matrixtyp mit selbstdefinierten Slicearray-Varianten entsprechend.)

Geben Sie die Eigenwerte und Eigenvektoren möglichst nicht in der Funktion, sondern im Hauptprogramm aus.

Rechnen Sie folgende Beispiele:

(a) A = ( \begin{matrix}5 & 4 & 1 & 1 \\ 4 & 5 & 1 & 1 \\ 1 & 1 & 4 & 2 \\ 1 & 1 & 2 & 4 \end{matrix} ) (b) A = ( \begin{matrix}5 & 1 & −2 & 0 & −2 & 5 \\ 1 & 6 & −3 & 2 & 0 & 6 \\ −2 & −3 & 8 & −5 & −6 & 0 \\ 0 & 2 & −5 & 5 & 1 & −2 \\ −2 & 0 & −6 & 1 & 6 & −3 \\ 5 & 6 & 0 & −2 & −3 & 8 \end{matrix} )

Bearbeitungszeitraum: bis Donnerstag, 26.1.2023, 16^{00}