Lineare Gleichungssysteme haben jeweils die Form $A \cdot x = b$ wobei $A$ und $b$ gegeben und $x$ gesucht ist:
$$A = \left(
\begin{matrix}
a_{11} & a_{12} & \cdots & a_{1n} \\
a_{21} & a_{22} & \cdots & a_{2n} \\
\vdots & \vdots && \vdots \\
a_{n1} & a_{n2} & \cdots & a_{nn}
\end{matrix}
\right),
x = \left(
\begin{matrix}
x_1 \\
x_2 \\
\vdots \\
x_n
\end{matrix}
\right),
b = \left(
\begin{matrix}
b_1 \\
b_2 \\
\vdots \\
b_n
\end{matrix}
\right)$$
</div>
### Eigenschaften
- Gleich viele gesuchte Variablen $x_n$ wie Gleichungen $n$. Folglich:
- Die Matrix $A$ ist eine quadratische Matrix mit Dimensionen $n \times n$
- $A$ ist invertierbar
- $A$ hat eine Determinante $\det(A)$
### Dreiecks-Matrizen
<divclass="letters">
***$L$: Untere Dreiecksmatrix***
Eine Matrix, die in der oberen rechten Ecke nur den Wert $0$ und auf der Diagonale nur den Wert $1$ hat. Eine Untere Dreiecksmatrix hat also folgende Form:
$$L = \left(
\begin{matrix}
1 & 0 & 0 & \cdots & 0 \\
l_{21} & 1 & 0 & \cdots & 0 \\
l_{31} & l_{32} & 1 & \ddots & 0 \\
\vdots & \vdots & \ddots & \ddots & 0 \\
l_{n1} & l_{n2} & \cdots & l_{nn - 1} & 1
\end{matrix}
\right)$$
***$R$: Obere Dreiecksmatrix***
Eine Matrix, die unten links von der Diagonale nur den Wert $0$ beinhaltet. Eine Obere Dreiecksmatrix hat dementsprechend folgende Form:
$$R = \left(
\begin{matrix}
r_{11} & r_{12} & r_{13} & \cdots & r_{1n} \\
0 & r_{22} & r_{23} & \cdots & r_{2n} \\
0 & 0 & r_{33} & \cdots & r_{3n} \\
\vdots & \vdots & \ddots & \ddots & \vdots \\
0 & 0 & \cdots & 0 & r_{nn}
\end{matrix}
\right)$$
</div>
***Code-Beispiele:***
_Umwandlung in $R$-Matrix:_
```py
for i in range(n):
if A[i, i] == 0:
index = -1
for j in range(i + 1, n):
if A[j, i] > 0:
index = j
if index == -1:
raise Exception("Invalid Matrix")
else:
# Swap lines
A[[i, index]] = A[[index, i]]
for j in range(i + 1, n):
factor = A[j, i] / A[i, i]
A[j] = A[j] - (factor * A[i])
```
### Der Gauss-Algorithmus
Der Gauss-Algorithmus basiert darauf, dass ein lineares Gleichungssystem leicht lösbar ist, falls $A$ eine obere Dreiecksmatrix ist. $A$ muss also hierfür die Form einer oberen Dreiecksmatrix $R$ haben.
<divclass="formula">
***Gauss-Algorithmus:***
$$x_i = \frac{b_i - \sum_{j = i + 1}^n{a_{ij} \cdot x_j}}{a_{ii}}, i = n, n - 1, \dots, 1$$
</div>
Um den Gauss-Algorithmus anzuwenden, muss die Matrix $A$ erst in eine $R$-Matrix umgewandelt werden. Dies funktioniert wie folgt:
1. Mit $i$ von $1$ bis $n$
2. Falls $a_{ii}$ den Wert $0$ hat:
1. Mit $j$ von $i + 1$ bis $n$
2. Prüfen, ob $a_{ji}$ einen höheren Wert als $0$ hat
- Falls Zeile gefunden wurde:
- $a_{i}$ mit $a_{j}$ tauschen
- $b_{i}$ mit $b_{j}$ tauschen
- Sonst beenden: ungültige Matrix
3. Mit $j$ von $i + 1$ bis $n$
1. $a_k = a_k - \frac{a_{ki}}{a_{ii}} \cdot a_i$
2. $b_k = b_k - \frac{a_{ki}}{a_{ii}} \cdot b_i$
***Code-Beispiel:***
```py
from numpy import array, zeros
def gaussMethod(A, b):
A = array(A)
n = A.shape[0]
A = A.reshape((n, n))
b = array(b).reshape((n))
result = zeros(n)
# Convert to R-Matrix
for i in range(n):
maxIndex = i
for j in range(i + 1, n):
if A[j, i] > A[maxIndex, i]:
maxIndex = j
# Swap lines
A[[i, maxIndex]] = A[[maxIndex, i]]
b[[i, maxIndex]] = b[[maxIndex, i]]
for j in range(i + 1, n):
factor = A[j, i] / A[i, i]
A[j] = A[j] - (factor * A[i])
b[j] = b[j] - (factor * b[i])
# Calculate result
for index in range(n, 0, -1):
i = index - 1
value = b[i]
for j in range(i, n):
value = value - A[i, j] * result[j]
result[i] = value / A[i, i]
return result.reshape((n, 1))
```
### Fehlerfortpflanzung und Pivotisierung
- Da beim Umwandeln einer Matrix $A$ in die $R$-Form Zeilen in jedem Schritt mit dem Faktor $\lambda = \frac{a_{ji}}{a_{ii}}$ multipliziert werden, vergrössert sich der Schritt immer um $|\lambda|$
- $\lambda$ kann klein gehalten werden, indem Zeilen der Grösse nach sortiert werden
- In den Code-Beispielen ist dies bereits berücksichtigt
### Determinanten-Bestimmung
Die Determinante einer Matrix $A$ lässt sich einfach berechnen, sobald sie in die $R$-Form gebracht wurde mit folgender Formel:
<divclass="formula">
Determinanten-Bestimmung mit Matrix $\tilde{A}$ (die Matrix $A$ in der $R$-Form):
$$\det(A) =
(-1)^l \cdot \det(\tilde{A}) =
(-1)^l \cdot \prod_{i = 1}^n{\tilde{a_{ii}}}$$
</div>
***Code-Beispiel:***
```py
from numpy import array
def det(A):
l = 0
n = A.shape[0]
A = A.reshape((n, n))
# Convert to R-Matrix
for i in range(n):
maxIndex = i
for j in range(i + 1, n):
if A[j, i] > A[maxIndex, i]:
maxIndex = j
# Swap lines
A[[i, maxIndex]] = A[[maxIndex, i]]
l = l + 1
for j in range(i + 1, n):
factor = A[j, i] / A[i, i]
A[j] = A[j] - (factor * A[i])
result = 1
for i in range(n):
result = result * A[i, i]
return (-1 ** l) * result
```
### Die $LR$-Zerlegung
In der $LR$-Zerlegung wird die Matrix $A$ in die Matrizen $L$ und $R$ aufgeteilt, sodass $A = L \cdot R$ gilt.
Alternative Namen dieses Vorgangs sind ***$LR$-Faktorisierung*** und $LU$-decomposition.
<divclass="formula">
Für in $L$ und $R$ zerlegte Matrizen gilt:
$$A \cdot x = b$$
und
$$A \cdot x = L \cdot R \cdot x = L \cdot y = b$$
Aufwand: Berechnung der $LR$-Zerlegung mit Gauss-Algorithmus benötigt ca. $\frac{2}{3}n^3$ Punktoperationen.
</div>
Falls Zeilenvertauschungen stattfinden, entsteht bei der $LR$-Zerlegung eine zusätzliche Permutations-Matrix $P$.
<divclass="formula">
Für $L$ und $R$ zerlegte Matrizen mit Permutation $P$ gilt:
$$P \cdot A = L \cdot R$$
$$L \cdot y = P \cdot b$$
$$R \cdot x = y$$
</div>
Das Verfahren für die $LR$-Zerlegung ist identisch zu den Schritten bei der Umwandlung in eine $R$-Matrix. Jedoch wird jeweils der Wert $l_{ji}$ in der (zu Beginn) leeren Matrix $L$ mit dem im aktuellen Eliminationsschritt gesetzt. Zudem muss bei Vertauschungen die Permutations-Matrix $P$ entsprechend angepasst werden:
***Code-Beispiel:***
```py
from numpy import array, identity, zeros
def decomposite(A):
l = 0
n = A.shape[0]
R = A.reshape((n, n))
L = zeros((n, n))
P = identity((n, n))
# Convert to LR-Matrix
for i in range(n):
maxIndex = i
for j in range(i + 1, n):
if A[j, i] > A[maxIndex, i]:
maxIndex = j
# Swap lines
Pn = identity((n, n))
A[[i, maxIndex]] = A[[maxIndex, i]]
Pn[[i, maxIndex]] = Pn[[maxIndex, i]]
P = P * Pn
for j in range(i + 1, n):
factor = R[j, i] / R[i, i]
L[j, i] = factor
R[j] = R[j] - (factor * R[i])
result = 1
for i in range(n):
result = result * R[i, i]
return [L, R, P]
```
Wenn die $LR$-Zerlegung, wie in diesem Code, Zeilenaustausch und das Berechnen von $P$ involviert, spricht man von einer $LR$-Zerlegung mit **Spaltenmaximum-Strategie**.
***Vorgang:***
1. Gemäss vorhergehender Beschreibung und Code-Beispiel die Matrizen $L$ und $R$ berechnen
2. Mit Hilfe des Gauss-Algorithmus $L \cdot y = P \cdot b$ nach $y$ auflösen
3. Mit Hilfe des Gauss-Algorithmus $R \cdot x = y$ nach $x$ auflösen
Beispiel der Berechnung einer Housholder-Matrix zur ersten Spalte der Matrix $A$.
> Für die Berechnung wird ein Einheitsvektor $e$ benötigt, welcher genauso viele Werte hat, wie die Matrix Dimensionen. Ein Einheitsvektor hat im ersten Feld den Wert $1$ und in allen anderen Feldern der Wert $0$.
>
> Für eine Matrix $A$ mit der Dimension $n = 3$ lautet der Einheitsvektor $e$ also wie folgt:
> $$e = \left(\begin{matrix}
> 1 \\
> 0 \\
> 0
> \end{matrix}\right)
1. Vektor $v$ bestimmen
$$v = a_1 + sign(a_{11}) \cdot |a_1| \cdot e$$
2. Vektor normieren:
$$u = \frac{1}{|v|} \cdot v =
\frac{1}{\sqrt{1^2 + 2^2 + 3^2}} \cdot
\left(\begin{matrix}
1 \\
2 \\
3
\end{matrix}\right) =
\frac{1}{\sqrt{14}} \cdot
\left(\begin{matrix}
1 \\
2 \\
3
\end{matrix}\right)$$
2. Die Housholder-Matrix $H = I_n - 2 \cdot u \cdot u^T$ berechnen.
$$H =
\left(\begin{matrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
\end{matrix}\right) -
2 \cdot \frac{1}{\sqrt{14}} \cdot
\left(\begin{matrix}
1 \\
2 \\
3
\end{matrix}\right) \cdot
\frac{1}{\sqrt{14}} \cdot
\left(\begin{matrix}
1 & 2 & 3
\end{matrix}\right) \\
H =
\left(\begin{matrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
\end{matrix}\right) -
2 \cdot \frac{1}{14} \cdot
\left(\begin{matrix}
1 & 2 & 3 \\
2 & 4 & 6 \\
3 & 6 & 9
\end{matrix}\right) =
-\frac{1}{7} \cdot
\left(\begin{matrix}
-6 & 2 & 3 \\
2 & -3 & 6 \\
3 & 6 & 2
\end{matrix}\right)$$
<divclass="letters">
- $H$: Housholder-Matrix
- $I$: Identitäts-Matrix
- $n$: Anzahl Dimensionen der Matrix
</div>
#### Vorgang
Im Rahmen des Vorgangs entspricht $A_1$ der Matrix $A$.
Die $QR$-Zerlegung kann folgendermassen durchgeführt werden:
Die Konditionszahl $cond(A)$ einer Matrix $A$ berechnet sich wie folgt:
$$cond(A) = ||A|| \cdot ||A^{-1}||$$
Eine hohe Konditionszahl $cond(A)$ bedeutet, dass kleine Fehler im Vektor $b$ zu grossen Fehlern im Ergebnis $x$ führen können. In diesem Fall ist eine Matrix schlecht konditioniert.
</div>
<divclass="formula">
***Fehlerabschätzung von Matrizen mit Fehlern:***
Sollte auch die Matrix $A$ fehlerhaft sein (fehlerhafte Matrix $\tilde{A}$), gilt der nachstehende Satz unter folgender Bedingung:
Die Ordnung $O(n)$ der zuvor genannten Verfahren entspricht der höchsten Potenz von $n$, welche in der Formel zur Berechnung des Aufwands vorkommt.
Das bedeutet also folgendes:
Ordnung von Gauss-Elimination, $LR$-Zerlegung und $QR$-Zerlegung:
$$O(n^3)$$
</div>
## Iterative Verfahren zur Lösung von Gleichungssystemen
### $LDR$-Zerlegung
Für die $LDR$-Zerlegung wird die Matrix $A$ in drei Matrizen $L$, $D$ und $R$ aufgeteilt, wobei $L$ eine untere Dreiecksmatrix, $D$ eine Diagonalmatrix und $R$ eine obere Dreiecksmatrix ist. Das bedeutet:
$$A = L + D + R$$
Mit
$$L = \left(
\begin{matrix}
0 & 0 & 0 & \cdots & 0 \\
a_{21} & 0 & 0 & \cdots & 0 \\
a_{31} & a_{32} & 0 & \cdots & 0 \\
\vdots & \vdots & \ddots & \ddots & \vdots \\
a_{n1} & a_{n2} & \cdots & a_{nn - 1} & 0
\end{matrix}
\right)$$
$$D = \left(
\begin{matrix}
a_{11} & 0 & 0 & \cdots & 0 \\
0 & a_{22} & 0 & \cdots & 0 \\
0 & 0 & a_{33} & \cdots & 0 \\
\vdots & \vdots & \ddots & \ddots & \vdots \\
0 & 0 & \cdots & 0 & a_{nn}
\end{matrix}
\right)$$
$$R = \left(
\begin{matrix}
0 & a_{12} & a_{13} & \cdots & a_{1n} \\
0 & 0 & a_{23} & \cdots & a_{2n} \\
0 & 0 & 0 & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & a_{n-1,n} \\
0 & 0 & \cdots & 0 & 0
\end{matrix}
\right)$$
> ***Wichtig:***
> Hierbei handelt es sich nicht um $L$ und $R$ aus der $LR$-Zerlegung!
### Jacobi-Verfahren
Das Jacobi-Verfahren ist ein iteratives Verfahren, welches nach jeder Iteration näher mit der tatsächlichen Lösung $x$ konvergiert.
Das Jacobi-Verfahren ist auch bekannt als **Gesamtschrittverfahren**.
<divclass="formula">
***Jacobi-Verfahren:***
Zunächst beginnt man mit $x^{(0)}$ als ein Vektor, der nur aus $0$en besteht.
Das Gauss-Seidel-Verfahren konvergiert schneller als das Jacobi-Verfahren.
Da für die Berechnung des Jacobi-Verfahrens für die Berechnung von $x_2$ auch Werte von $x_1$ verwendet werden, können die Werte direkt aus der aktuellen Iteration $k$ wiederverwendet werden, um den Vorgang schneller konvergieren zu lassen.
Das Gauss-Seidel-Verfahren wird auch **Einzelschrittverfahren** genannt.
$$x^{(n + 1)} = B \cdot x^{(n)} + c =: F(x^{(n)})$$
Beispiele für solche Fixpunkt-Iterationen sind das Jacobi- oder das Gauss-Seidel-Verfahren.
Falls folgendes gegeben ist: $\tilde{x} = B \cdot \tilde{x} + c = F(\tilde{x})$, dann gilt:
- $\tilde{x}$ ist ein anziehender Fixpunkt, falls $||B|| <1$
- $\tilde{x}$ ist ein abstossender Fixpunkt, falls $||B|| > 1$
</div>
<divclass="formula">
***Abschätzungen:***
Gegeben sei eine Fixpunkt-Iteration:
$$x^{(n + 1)} = B \cdot x^{(n)} + c =: F(x^{(n)})$$
Für Fixpunkt-Iterationen bei denen $x^{(k)}$ gegen $\tilde{x}$ konvergiert (gemäss oben stehender Formel "Anziehung/Abstossung"), gelten folgende Abschätzungen:
besitzt in der Menge $\mathbb{C}$ der komplexen Zahlen genau $n$ Lösungen.
</div>
<divclass="formula">
***Ziehen der Wurzel einer komplexen Zahl:***
Die Gleichung für das Ziehen einer Wurzel $n$ der komplexen Zahl $a$ lautet: $z^n = a$.
Für die Lösung dieser Gleichung existieren genau $n$ verschiedene Lösungen in der Menge $\mathbb{C}$:
$$z_k = r \cdot (\cos(\varphi_k + i \cdot \sin(\varphi_k)) = r \dot e^{i \cdot \varphi_k}$$
für $k = 0, 1, 2, \dots, n - 1$:
mit
$$r = \sqrt[n]{r_0}$$
$$\varphi_k = \frac{\varphi + k \cdot 2 \cdot \pi}{n}$$
Die Bildpunkte der Ergebnisse liegen in der komplexen Zahlenebene auf einem Kreis um den Nullpunkt mit dem Radius $r = \sqrt[n]{r_0}$ und bilden die Ecken eines regelmässigen $n$-Ecks.
</div>
Visualisierung des Ziehens der $6$-ten Wurzel einer komplexen Zahl $z$:
Die Vektor-Iteration, auch **von-Mises-Iteration** genannt, erlaubt das Bestimmen des grössten Eigenwertes $\lambda$ einer diagonalisierbaren Matrix $A$.
<divclass="formula">
***Spektral-Radius:***
Der Spektral-Radius $\rho(A)$ definiert den höchsten Eigenwert der Matrix $A$:
$$\rho(A) = \max\{|\lambda|\; | \; \lambda \text{ ist ein Eigenwert von }A \in \mathbb{R}^{n \times n}\}$$
</div>
<divclass="formula">
Sei $A$ eine diagonalisierbare Matrix mit den Eigenwerten $\lambda_1, \dots, \lambda_n$ wobei $\lambda_1$ betragsmässig am höchsten ist:
