Commit c7934f55 authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Bessere Beschreibung iterativer Löser

parent 1489c695
Pipeline #6192 passed with stage
in 3 minutes and 4 seconds
......@@ -43,3 +43,10 @@ year = {2013}
volume = {49},
pages = {409--436}
}
@Book{hackbusch:1994,
title = {Iterative Solution of Large Sparse Systems of Equations},
publisher = {Springer},
year = {1994},
author = {Wolfgang Hackbusch}
}
......@@ -552,27 +552,22 @@ Direkte Verfahren wie Gauß-Elimination oder Cholesky-Zerlegung sind i.A. zu teu
Iterative Verfahren konstruieren eine Folge $\{ x^k \}_{k \in \N}$,
die (hoffentlich) gegen die Lösung $x^*$ von $Ax=b$ konvergiert.
\medskip
\bigskip
Bei linearen iterativen Verfahren ist die Abbildung
Für \emph{lineare} iterative Verfahren~\cite{dahmen_reusken:2008} schreibt man
$Ax=b$ als Fixpunktgleichung
\begin{equation*}
e^k \colonequals x^k - x^*
\qquad \mapsto \qquad
e^{k+1} \colonequals x^{k+1} - x^*
x=x+C(b-Ax)
\end{equation*}
linear.
\bigskip
mit einer nichtsingulären Matrix $C \in \R^{n \times n}$.
Der normale Ansatz ist~\cite{dahmen_reusken:2008}:
\medskip
Schreibe $Ax=b$ als Fixpunktgleichung
\begin{equation*}
x=x+C(b-Ax)
\end{equation*}
mit $C \in \R^{n \times n}$ nichtsingulär. Die Größe $r(x) \colonequals b - Ax$
Die Größe $r(x) \colonequals b - Ax$
nennt man \emph{Residuum}. Die Matrix $C$ heißt \emph{Vorkonditionierer}.
\bigskip
Dafür machen wir jetzt eine Fixpunktiteration:
\begin{equation*}
x^{k+1}\colonequals x^k+C (b-Ax^k)=(I-CA)x^k+Cb,
......@@ -580,44 +575,29 @@ Dafür machen wir jetzt eine Fixpunktiteration:
k = 0, 1, 2, \dots
\end{equation*}
Unter welchen Umständen konvergiert dieses Verfahren?
\medskip
Der Fehler im $k$-ten Iterationsschritt ist $e^k=x^k-x^*$. Es gilt
Der Fehler im $k$-ten Iterationsschritt ist $e^k \colonequals x^k-x^*$. Es gilt
\begin{equation*}
e^{k+1}=x^{k+1}-x^*
=
\underbrace{(I-CA)}_{\text{Iterationsmatrix}}e^k.
\end{equation*}
Also gilt
\begin{equation}
\label{equa:eqerrlin}
e^k=(I-CA)^ke^0
\qquad
\forall k=0,1,2,\ldots
\end{equation}
Das erklärt den Namen: Die Abbildung $e^k \mapsto e^{k+1}$ ist linear.
\begin{definition}
Die Matrix $I - CA$ heißt \emph{Iterationsmatrix} der Methode.
\end{definition}
\subsubsection{Konvergenz:}
Die Fehlerfortpflanzung \eqref{equa:eqerrlin} ist tatsächlich linear.
\bigskip
Für die Konvergenz gilt der folgende wichtige Satz.
Für jede Vektornorm mit assoziierter Matrixnorm erhält man
\begin{theorem}
Sei $\norm{\cdot}$ eine submultiplikative Matrixnorm.
Das Verfahren konvergiert für jeden Startwert $x^0 \in \R^n$ gegen die Lösung
von $Ax = b$, wenn $\norm{I - CA} < 1$.
\end{theorem}
Da $\rho (B) \leq \norm{B}$ für jede submultiplikative Matrixnorm gilt sogar
[Denn: Sei $v$ ein Eigenvektor von $B$ zum Eigenwert $\lambda$. Dann ist
Sei $\rho$ der Spektralradius einer Matrix. Für jede submultiplikative Matrixnorm
gilt $\rho (B) \leq \norm{B}$, denn für jeden Eigenvektor $v$ von $B$
zum Eigenwert $\lambda$ gilt
\begin{equation*}
\abs{\lambda}\norm{v}
=
......@@ -627,13 +607,22 @@ Da $\rho (B) \leq \norm{B}$ für jede submultiplikative Matrixnorm gilt sogar
\le
\norm{B}\norm{v}.
\end{equation*}
Deshalb gilt $\abs{\lambda} \le \norm{B}$ für alle Eigenwerte $\lambda$ von $B$.]
Damit erhält man:
\begin{theorem}
Sei $\rho (I-CA)$ der Spektralradius von $I-CA$. Das Verfahren konvergiert für jeden Startwert $x^0 \in \R$ gegen die Lösung von $Ax=b$ genau dann, wenn $\rho(I-CA)<1$.
\end{theorem}
Je nach Wahl des Vorkonditionierers $C$ erhält man unterschiedliche Verfahren:
Je nach Wahl des Vorkonditionierers $C$ erhält man unterschiedliche Verfahren.
Es ergibt sich folgendes Dilemma:
\begin{enumerate}
\item $C$ soll $A^{-1}$ möglichst gut approximieren,
\item die Operation $y \mapsto Cy$ soll möglichst billig sein.
\end{enumerate}
Seien $D$, $L$, $R$ die Diagonale bzw.\ der linke und rechte Dreiecksteil von $A$.
Man erhält die folgenden Standardverfahren:
\begin{itemize}
\item Das Jacobi-Verfahren
\begin{equation*}
......@@ -647,15 +636,10 @@ Je nach Wahl des Vorkonditionierers $C$ erhält man unterschiedliche Verfahren:
\item Das symmetrische Gauß--Seidel-Verfahren
\begin{equation*}
C = (D-L)^{-1}+(D-R)^{-1}-(D-R)^{-1}A(D-L)^{-1}
C = (D+L)^{-1}+(D+R)^{-1}-(D+R)^{-1}A(D+L)^{-1}
\end{equation*}
\end{itemize}
Es ergibt sich folgendes Dilemma:
\begin{enumerate}
\item $C$ soll $A^{-1}$ möglichst gut approximieren
\item Die Operation $y \mapsto Cy$ soll möglichst billig sein
\end{enumerate}
\emph{Beachte:} Die Matrix $C$ wird nie explizit ausgerechnet!
......@@ -664,35 +648,35 @@ Es ergibt sich folgendes Dilemma:
Wir betrachten dieses Verfahren etwas genauer. Man kann es auch
\glqq direkt\grqq motivieren.
Wir nehmen im Folgenden an, dass $a_{ii} \neq 0$ für alle $i=1,\ldots,n$.
Wir nehmen im Folgenden an, dass $A_{ii} \neq 0$ für alle $i=1,\ldots,n$.
\medskip
Betrachte das lineare Gleichungssystem:
\begin{align*}
a_{11} x_1+a_{12} x_2+\ldots +a_{1m} x_m & = b_1 \\
a_{12} x_1+a_{22} x_2+\ldots +a_{2m} x_m & = b_2 \\
& \vdots
A_{11} x_1+A_{12} x_2+\ldots +A_{1m} x_m & = b_1 \\
A_{12} x_1+A_{22} x_2+\ldots +A_{2m} x_m & = b_2 \\
& \;\; \vdots
\end{align*}
\emph{Idee:} Löse $i$-te Zeile nach $x_i$ auf für $i=1,\ldots,n$.
\emph{Idee:} Löse $i$-te Zeile nach $x_i$ auf für $i=1,\ldots,n$:
\begin{align*}
x_1 & =\frac{1}{a_{11}} \left(b_1-a_{12} x_2-\ldots -a_{1n} x_n \right) \\
x_2 & =\frac{1}{a_{22}} \left(b_2-a_{21} x_1-\ldots -a_{2n} x_n \right) \\
& \vdots
x_1 & =\frac{1}{A_{11}} (b_1-A_{12} x_2-\ldots -A_{1n} x_n ) \\
x_2 & =\frac{1}{A_{22}} (b_2-A_{21} x_1-\ldots -A_{2n} x_n ) \\
& \;\; \vdots
\end{align*}
Mache daraus ein iteratives Verfahren:
\begin{align*}
x_1^{k+1} & =\frac{1}{a_{11}} \left(b_1-a_{12}x_2^k-\ldots -a_{1n}x_n^k \right) \\
x_2^{k+1} & =\frac{1}{a_{22}} \left(b_2-a_{21}x_1^k-\ldots -a_{2n}x_n^k \right) \\
& \vdots
x_1^{k+1} & =\frac{1}{A_{11}} (b_1-A_{12}x_2^k-\ldots -A_{1n}x_n^k) \\
x_2^{k+1} & =\frac{1}{A_{22}} (b_2-A_{21}x_1^k-\ldots -A_{2n}x_n^k) \\
& \;\; \vdots
\end{align*}
Für $i=1,\ldots,n$
\begin{equation*}
x_i^{k+1}
=
\frac{1}{a_{ii}} \bigg(b_i-\sum_{\substack{j=1\\j\neq i}}^n a_{ij}x_i^k \bigg)
\frac{1}{A_{ii}} \bigg(b_i-\sum_{\substack{j=1\\j\neq i}}^n A_{ij}x_i^k \bigg)
=
x^k_i + \frac{1}{a_{ii}} \bigg(b_i-\sum_{j=1}^n a_{ij}x_i^k \bigg).
x^k_i + \frac{1}{A_{ii}} \bigg(b_i-\sum_{j=1}^n A_{ij}x_i^k \bigg).
\end{equation*}
Beachte:
\begin{itemize}
......@@ -706,17 +690,29 @@ Beachte:
$C = D^{-1}$ ist.
\end{exercise}
\todo[inline]{Hier muss direkt das gedämpfte Verfahren eingeführt werden.}
Häufig erweitert man das Jacobi-Verfahren zum \emph{gedämpften} Jacobi-Verfahren.
Dabei wählt man einen Parameter $\eta > 0$ und definiert
\begin{equation*}
x_i^{k+1}
=
x^k_i + \frac{\eta}{A_{ii}} \bigg(b_i-\sum_{j=1}^n A_{ij}x_i^k \bigg)
\qquad \text{bzw.} \qquad
x^{k+1} = x^k + \eta D^{-1} ( b - Ax^k).
\end{equation*}
\subsubsection{Konvergenz}
Wir wenden das bekannte Konvergenzkriterium an:
\begin{theorem}
Das Jacobi Verfahren konvergiert genau dann, wenn $\rho \left(I-D^{-1}A \right) <1$.
Das Jacobi Verfahren konvergiert genau dann, wenn $\rho \left(I- \eta D^{-1}A \right) <1$.
\end{theorem}
Leider gilt diese Bedingung nicht immer.
Ein paar direkte Charakterisierungen findet man z.B.\ bei \citet{dahmen_reusken:2008}.
Das Verfahren konvergiert also immer, wenn man nur $\eta$ klein genug wählt.
Allerdings wird das Verfahren für kleine $\eta$ auch immer langsamer.
\medskip
Ein paar direkte Charakterisierungen für den Fall $\eta = 1$ findet man z.B.\ bei \citet{dahmen_reusken:2008}.
Wir greifen die Konvergenztheorie des Jacobi-Verfahrens in Kapitel~\ref{sec:jacobi_smoother}
wieder auf.
......@@ -741,12 +737,12 @@ mit Dreiecks- oder Viereckselementen. Eine Zeile von $A$ ist dann
Jacobi-Verfahren: $D=4h^{-2}I$.
Wir versuchen jetzt, den Spektralradius von $I-CA = I - D^{-1}A$ abzuschätzen:
Wir versuchen jetzt, den Spektralradius von $I-CA = I - \eta D^{-1}A$ abzuschätzen:
\begin{align*}
\rho (I-D^{-1}A)
& = \sup_{x \neq 0} \abs[\bigg]{\frac{x^T \left(I-D^{-1}A \right) x}{\norm{x}^2}} \\
& = \sup_{x \neq 0} \abs[\bigg]{1-\frac{x^T D^{-1}Ax}{\norm{x}^2}} \\
& = \sup_{x \neq 0} \abs[\bigg]{1-\frac{1}{4}h^2 \frac{x^TAx}{\norm{x}^2}}.
\rho (I - \eta D^{-1}A)
& = \sup_{x \neq 0} \abs[\bigg]{\frac{x^T (I- \eta D^{-1}A) x}{\norm{x}^2}} \\
& = \sup_{x \neq 0} \abs[\bigg]{1- \eta \frac{x^T D^{-1}Ax}{\norm{x}^2}} \\
& = \sup_{x \neq 0} \abs[\bigg]{1-\frac{1}{4}\eta h^2 \frac{x^TAx}{\norm{x}^2}}.
\end{align*}
Der Bruch ist gerade der \emph{Rayleigh-Quotient}
\begin{equation*}
......@@ -755,25 +751,25 @@ Der Bruch ist gerade der \emph{Rayleigh-Quotient}
Für symmetrische $A$ gilt $\lambda_\text{min}(A) \le R(A,x) \le \lambda_\text{max}(A)$, und diese Schranken werden
angenommen, wenn $x$ entsprechende Eigenvektoren sind.
\medskip
Daraus folgt dass
\begin{equation*}
\rho (I-D^{-1}A )
\begin{align*}
\rho (I - \eta D^{-1}A )
& =
\sup \left\lbrace \abs[\Big]{1- \frac{1}{4} \eta h^2 \lambda} \; : \; \text{$\lambda$ Eigenwert von $A$} \right\rbrace \\
%
& =
1-2 \eta \sin^2 \left( \frac{1}{2} \pi h \right).
\end{align*}
Mit $\eta = 1$ erhält man:
\todo[inline]{Zeige den allgemeinen Fall!}
\begin{align*}
\rho (I - \eta D^{-1}A )
=
\sup \left\lbrace \abs[\Big]{1- \frac{1}{4} h^2 \lambda} \; : \; \text{$\lambda$ Eigenwert von $A$} \right\rbrace.
\end{equation*}
Für die spezielle Matrix können wir die Eigenwerte ausrechnen. Diese sind
alle nichtnegativ.
Der größte Eigenwert von $A$ ist~\cite[Kapitel~12.3.3]{dahmen_reusken:2008}:
\cos (\pi h).
\end{align*}
\begin{equation*}
\lambda =\frac{8}{h^2} \sin^2 \left( \frac{1}{2} \pi h \right)
\end{equation*}
Es folgt, dass
\begin{equation*}
\rho (I-D^{-1}A)=1-2 \sin^2 \left( \frac{1}{2} \pi h \right)=\cos (\pi h)
\end{equation*}
Taylor-Entwicklung:
\begin{equation*}
\cos \left( \pi h \right)=1-\frac{1}{2} \pi^2 h^2+\ldots
......@@ -803,9 +799,11 @@ Da $\frac{\norm{e^k}}{\norm{e^0}} \approx \rho^k$ erhält man
\end{equation*}
da $\ln x \approx (x-1)-\frac{1}{2} (x-1)^2 + \ldots$ ist.
Für ein doppelt so feines Gitter braucht man viermal so viele Iterationen
Für ein doppelt so feines Gitter braucht man also viermal so viele Iterationen
(und diese sind natürlich auch noch teurer.).
\medskip
Das Jacobi-Verfahren scheint also kein sonderlich gutes Verfahren zu sein.
\subsection{Das Gauß-Seidel-Verfahren}
......@@ -840,9 +838,10 @@ Wir erwarten bessere Konvergenzeigenschaften als für das Jacobi-Verfahren.
Und in der Tat:
\begin{theorem}
Das Gauß-Seidel-Verfahren konvergiert, wenn \begin{itemize}
\item $A$ symmetrisch und positiv definit ist und/oder
\item $A$ irreduzibel und diagonal dominant ist.
Das Gauß-Seidel-Verfahren konvergiert, wenn
\begin{itemize}
\item $A$ symmetrisch und positiv definit ist und/oder
\item $A$ irreduzibel und diagonaldominant ist.
\end{itemize}
\end{theorem}
\subsubsection{Konvergenzgeschwindigkeit}
......@@ -938,8 +937,8 @@ Sei $\kappa$ die Kondition der Matrix $A$. Dann gilt nach $k$ Schritten des CG-V
\item Dabei wird wieder die Gleichung $0 = b - Ax$ durch $0 = C(b-Ax)$ ersetzt.
\item $C$ muss so gewählt werden dass $\kappa(CA)$ klein ist, aber $Cv$, $v \in \R^n$
billig ausgerechnet werden kann.
\item $C$ muss so gewählt werden dass $\kappa(CA)$ klein ist, dass aber gleichzeitig
$Cv$ für $v \in \R^n$ billig ausgerechnet werden kann.
\item Die Konstruktion passender $C$ bildet ein großes Gebiet innerhalb
der Numerik.
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment