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

Programmieren I für Studierende der Mathematik

Aufgabe 7

Eine symmetrische, positiv definite Matrix A = (a_{ij})_{i,j=0,...,n−1} ∈ ℝ^{n×n} ist darstellbar als Produkt A = LL^{T}, wobei L = (l_{ij})_{i,j=0,...,n−1} eine untere Dreiecksmatrix mit positiven Hauptdiagonalelementen ist (Cholesky-Zerlegung). L kann spaltenweise mit folgenden Formeln berechnet werden.

\begin{matrix}l_{jj} = \sqrt{a_{jj} − ∑_{k=0}^{j−1} l_{jk}^{2}} \\ l_{ij} = \frac{1}{l_{jj}} (a_{ij} − ∑_{k=0}^{j−1} l_{ik}l_{jk}) i = j + 1, …, n − 1 \end{matrix} \} j = 0, …, n − 1.

Ein lineares Gleichungssystem Ax = b kann dann bei gegebenem b = (b_{i})_{i=0,...,n−1} ∈ ℝ^{n} und gesuchtem x = (x_{i})_{i=0,...,n−1} ∈ ℝ^{n} durch Lösen von Ly = b und L^{T} x = y aufgelöst werden, was zu den Formeln

\begin{matrix}y_{i} = \frac{1}{l_{ii}} (b_{i} − ∑_{k=0}^{i−1} l_{ik}y_{k}) & i = 0, …, n − 1 \\ x_{i} = \frac{1}{l_{ii}} (y_{i} − ∑_{k=i+1}^{n−1} l_{ki}x_{k}) & i = n − 1, …, 0. \end{matrix}

führt. Das Programm soll mit einer Fehlermeldung abgebrochen werden, wenn der Radikant bei

der Berechnung von l_{jj} kleiner oder gleich 0 oder kleiner als ε · Spur(A) ist (Spur(A) := ∑_{j=0}^{n−1} a_{jj}).

Erstellen Sie eine Funktion cholesky, die die Matrixgleichung Ax = b auf diese Weise löst und im Erfolgsfall den Vektor x als Ergebnis liefert. Der Funktion sollen lediglich die Größen A, b und ε als Argumente übergeben werden.

Das Einlesen von n, A und b (mit Kontrollausgabe!) und die Ausgabe von x soll im Hauptprogramm erfolgen.

Rechnen Sie folgende Beispiele (mit ε = 10^{−10})

(a) A = ( \begin{matrix}1 & 2 & 4 \\ 2 & 13 & 23 \\ 4 & 23 & 77 \end{matrix} ) , b = ( \begin{matrix}−3 \\ −21 \\ −73 \end{matrix} ) (b) A = ( \begin{matrix}1 & 0 & 3 & 0 \\ 0 & 4 & 0 & 0 \\ 3 & 0 & 25 & 20 \\ 0 & 0 & 20 & 61 \end{matrix} ) , b = ( \begin{matrix}10 \\ 8 \\ 158 \\ 304 \end{matrix} )

(c) A = ( \begin{matrix}10 & 1 & 2 & 3 & 4 \\ 1 & 9 & −1 & 2 & −3 \\ 2 & −1 & 7 & 3 & −5 \\ 3 & 2 & 3 & 12 & 1 \\ 4 & −3 & −5 & 1 & 15 \end{matrix} ) , b = ( \begin{matrix}12 \\ −8 \\ 1 \\ 3 \\ 4 \end{matrix} )

Bearbeitungszeitraum: bis Freitag, 23.6.2023, 12^{00}