Verbessere die anschauliche Einführung des Mehrgitterverfahrens

skript-mehrgitter-sander.tex

@@ -930,154 +930,301 @@ Sei $\kappa$ die Kondition der Matrix $A$. Dann gilt nach $k$ Schritten des CG-V
  \item Eine Möglichkeit: Mehrgitter!
 \end{itemize}
 
+\section{Mehrgitterverfahren: Der anschauliche Zugang}
 
+\subsection{Glättung}
 
+Gauß--Seidel-Verfahren und Jacobi-Verfahren sind für typische Finite-Elemente-Matrizen
+extrem langsam.
 
-\subsubsection{Unvollständige Cholesky-Zerlegung (ICH,ILU,\ldots)}
+\medskip
+
+Sie haben allerdings eine wichtige Eigenschaft, die wir uns
+für Mehrgitterverfahren zunutze machen:
+\begin{itemize}
+ \item Die zwei Verfahren \emph{glätten}.
+\end{itemize}
 
+\emph{Beispiel:} Wir betrachten das Poisson-Problem
 
-Für die Wahl von $W$ sind extrem viele verschiedene Möglichkeiten vorgeschlagen worden.
+\begin{align*}
+ -\Delta u & = 1  \qquad \text{auf $\Omega$}, \\
+         u & = 0  \qquad \text{auf $\partial \Omega$,}
+\end{align*}
+auf dem Einheitsquadrat $\Omega = (0,1)^2$.
 
 \medskip
 
-Wir zeigen zwei wichtige Ansätze.
+Hier ist die Finite-Elemente-Lösung auf einem strukturierten Dreiecksgitter
+mit $16 \times 16 \times 2$ Elementen, dargestellt als Höhenplot:
+\begin{center}
+ \includegraphics[width=0.6\textwidth,trim=0 50 0 10]{gs_from_zero_iterate50}
+\end{center}
+
+Wir versuchen jetzt, das dazugehörige algebraische Problem $Ax = b$
+mit dem Gauß--Seidel-Verfahren zu lösen.
+
+\bigskip
+
+Eine Iteration ist:
+
+\begin{algorithm}[H]
+ \SetAlgoLined
+  % \caption{...}
+ \For{alle Zeilen $i$}
+ {
+  \begin{equation*}
+   x^{k+1}_i
+    \longleftarrow
+   \frac{1}{A_{ii}}
+     \bigg( b_i - \sum_{j=1}^{i-1} A_{ij} x^{k+1}_j - \sum_{j=1+1}^{n} A_{ij} x^k_j \bigg)
+  \end{equation*}
+ }
+\end{algorithm}
+
+Als Startiterierte $x^0$ wählen wir $x^0 = 0 \in \R^n$.
+
+\bigskip
+
+Hier sind einige der ersten Iterierten:
+
+\begin{center}
+ \begin{tabular}{cc}
+
+  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_zero_iterate0} &
+  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_zero_iterate1} \\
+
+  iterate $x^0$ & iterate $x^1$ \\
+
+  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_zero_iterate2} &
+  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_zero_iterate3} \\
+
+  iterate $x^2$ & iterate $x^3$ \\
+
+  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_zero_iterate10} &
+  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_zero_iterate50} \\
+
+  iterate $x^{10}$ & iterate $x^{50}$
+ \end{tabular}
+\end{center}
+
+Man erkennt dass das Verfahren gehen die Lösung konvergiert, aber es braucht
+tatsächlich recht lange.
+
+\bigskip
+
+Jetzt wollen wir zeigen was wir mit \glqq Gauß--Seidel glättet\grqq{} meinen.
 
 \medskip
 
-\emph{Idee:} Die beste Kondition bekäme man natürlich für $W=A$.
+Wir lösen nochmal die Gleichung, fangen jetzt aber von einem
+verrauschten Startwert an.
+
+\begin{center}
+ \begin{tabular}{cc}
+  \includegraphics[width=0.3\textwidth,trim=0 50 0 0]{gs_from_noise_iterate0} &
+  \includegraphics[width=0.3\textwidth,trim=0 50 0 0]{gs_from_noise_iterate1} \\
+
+  iterate $x^0$ & iterate $x^1$ \\
+
+  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_noise_iterate2} &
+  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_noise_iterate3} \\
+
+  iterate $x^2$ & iterate $x^3$ \\
+
+  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_noise_iterate10} &
+  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_noise_iterate50} \\
+
+  iterate $x^{10}$ & iterate $x^{50}$
+ \end{tabular}
+\end{center}
+
 \begin{itemize}
- \item Zur Lösung von $Wz=Az=r_k$ könnte man die Cholesky-Zerlegung $A=LDL^T$ berechnen.
+ \item Das Verfahren konvergiert immer noch.
 
- \item Das ist vermutlich keine gute Idee, denn wenn man die
- Cholesky-Zerlegung dann hat kann man direkt $Ax=b$ lösen,
- und braucht kein CG-Verfahren mehr.
+ \item Es scheint weder schneller noch langsamer geworden zu sein.
 
- \item ... aber ignorieren wir das mal...
+ \item Das Rauschen ist aber schon nach sehr wenigen Iterationen weitestgehend
+   verschwunden.
+\end{itemize}
 
- \item Im Zuge des vorkonditionierten CG-Verfahrens werden immer wieder
-  Gleichungssysteme mit der Matrix $W$ gelöst.  Es kann sich also
-  lohnen ein Zerlegung von $W$ anzufertigen, denn man kann sie
-  mehrfach anwenden.
 
- \item Selbst für dünnbesetzte $A$ ist $L$ aber vollbesetzt \newline
- $\implies$ Lösen mit Cholesky-Zerlegung ist zu teuer und braucht zu viel Speicher.
+Das klingt komisch:
+\begin{itemize}
+ \item Was würde passieren wenn die rechte Seite so gewählt wäre,
+   dass die Lösung des Poisson-Problems ein (festes) Rauschen wäre?
 
- \item Aber es würde ja reichen wenn $W \approx A$.
+ \item In so einem Fall könnte doch das Gauß--Seidel-Verfahren nicht glätten?
 \end{itemize}
 
-\begin{definition}
-Sei $A \in \R^{n \times n}$ symmetrisch und positiv definit. Eine unvollständige Cholesky-Zerlegung von $A$ ist
-$A \approx \tilde{L} \tilde{D} \tilde{L}^T$, wobei
+In Wirklichkeit ist es etwas anders:
+
 \begin{itemize}
- \item $\tilde{L}$ : normierte untere Dreiecksmatrix
- \item $\tilde{D}$ : Diagonalmatrix
- \item $\tilde{L}_{ij}=0$ falls $A_{ij}=0$
+ \item Das Gauß--Seidel-Verfahren glättet \emph{nicht} die Iterierten.
+
+ \item Das Gauß--Seidel-Verfahren glättet die \emph{Fehler}!
 \end{itemize}
-\end{definition}
 
-\bigskip
+\subsection{Fourier-Analyse}
+\label{sec:jacobi_smoother}
+
+\subsection{Zweigitterverfahren}
+
+Die zweite wichtige Einsicht ist:
 
-Vorkonditionierer: $W=\tilde{L} \tilde{D} \tilde{L}^T$
 \begin{itemize}
- \item häufig: $\kappa \left(W^{-1}A \right) \ll \kappa (A)$
- \item $\tilde{L}$ ist dünnbesetzt: $\implies \tilde{L} \tilde{D} \tilde{L}^Tz^k=r_k$ ist billig zu lösen.
+ \item Der Fehler $e^k \colonequals x^* - x^k$ löst eine lineare Gleichung
+  \begin{equation*}
+   A(x^k + e) = b \qquad \iff \qquad Ae = b - Ax^k \equalscolon r^k.
+  \end{equation*}
+
+  \item Da $e^k$ glatt ist reicht es, $Ae = r^k$ auf einem gröberen Gitter lösen.
 \end{itemize}
 
+Wir konstruieren uns deshalb ein zweites gröberes Gitter $\mathcal{T}_\text{grob}$ so,
+dass der dazugehörige Finite-Elemente-Raum $V_h^\text{grob}$ ein Teilraum
+von $V_h^\text{fein} = V_H$ ist.
 
-\begin{example}
-Poisson-Problem. CG vs.\ CG mit ILR-Vorkonditionierer
-\begin{equation*}
-\begin{tabular}{c|c|ccc}
-h & $\frac{1}{40}$ & $\frac{1}{80}$ & $\frac{1}{160}$ & $\frac{1}{320}$ \\
-\hline
-CG & 65 & 130 & 262 & 525 \\
-ILR-CG & 20 & 40 & 79 & 157
-\end{tabular}
-\end{equation*}
-Tabelle: Anzahl der Iterationen um das Residuum um den Faktor $1000$ zu reduzieren.
-\end{example}
+   \begin{center}
+    \begin{tikzpicture}[scale=0.35]
+     % Feines Gitter
+     \draw [step=0.5] (0,0) grid (8,8);
 
-Es sind viele Varianten möglich!
+     \foreach \i in {0.5, 1, ..., 8} {
+      \draw (0,\i) -- (\i,0);
+      \draw (\i,8) -- (8,\i);
+     }
 
-\subsubsection{Lineare Verfahren als Vorkonditionierer}
-\emph{Erinnerung:} Lineares Verfahren
-\begin{equation*}
- x^{k+1}=x^k+C (b-Ax^k)
-\end{equation*}
-Sei
-\begin{enumerate}
- \item $A$ symmetrisch und positiv definit
- \item $C$ so, dass $\rho (I-CA)<1$, d.h. das Verfahren konvergiert.
-\end{enumerate}
-Wir hatten $C$ als Approximation von $A^{-1}$ interpretiert.
+     % Grobes Gitter
+     \draw [step=1.0] (15,0) grid (23,8);
 
-\bigskip
+     \foreach \i in {1, 2, ..., 8} {
+      \draw (15,\i) -- (15+\i,0);
+      \draw (15+\i,8) -- (23,\i);
+     }
+    \end{tikzpicture}
 
-\emph{Idee}: Wähle $W^{-1}=C$ als Vorkonditionierer für CG.
+   \end{center}
 
-\medskip
+\begin{itemize}
+  \item Definiere Prolongationsoperator als die kanonische Injektion
+    \begin{equation*}
+     P : V_h^\text{grob} \to V_h^\text{fein}
+    \end{equation*}
 
-\emph{Praktische Umsetzung}: Für CG müssen Ausdrücke der Form $z = W^{-1} r_k = Cr_k$ berechnet werden.
+    \begin{center}
+     \includegraphics[width=0.3\textwidth]{coarse_function}
+     \hspace{0.1\textwidth}
+     \includegraphics[width=0.3\textwidth]{prolonged_coarse_function}
+    \end{center}
+\end{itemize}
+
+Mit $P$ kann man Funktionen vom groben auf das feine Gitter transferieren,
+    ohne einen Fehler zu machen.
 
 \medskip
 
-\emph{Problem}: $C$ ist i.\,A.\ nicht explizit gegeben.
+\begin{itemize}
+ \item Natürlich können wir keine Funktion vom feinen Gitter auf das grobe Gitter
+  übertragen, ohne einen Approximationsfehler zu machen.
 
-\smallskip
+ \item \emph{Aber:} Wir können Linearformen $\ell(\cdot) : V_h^\text{fein} \to \R$
+  auf den groben FE-Raum beschränken:
+  \begin{equation*}
+   \ell^\text{grob}(v_h) \colonequals \ell^\text{fein}(Pv_h)
+   \qquad
+   \forall v_h \in V_h^\text{grob}.
+  \end{equation*}
 
-\emph{Lösung:} Betrachte das lineare Gleichungssystem $Az=r_k$
-\begin{itemize}
- \item Setze $z^0=0 \in \R^n$
- \item Ein Schritt des linearen Verfahrens:
+ \item Das Residuum $r^k$ ist die algebraische Darstellung so einer Linearform:
+
+  Funktionenraum:
   \begin{equation*}
-   z^1=z^0+C (r_k-Az^0)=Cr_k
+   r^k(v_h) \colonequals \langle b,v_h \rangle - a(u_h^k,v_h)
+   \qquad
+   \forall v_h \in V_h^\text{fein}
   \end{equation*}
-\end{itemize}
 
-\bigskip
+  Algebraisch :
+  \begin{equation*}
+   (r^k)^T v \colonequals b^T v - (x^k)^T A v
+   \qquad
+   \forall v \in \R^n
+  \end{equation*}
 
-Viele lineare Verfahren können als Vorkonditionierer verwendet werden!
+ \item Die Einschränkung auf den groben Funktionenraum ist
+  \begin{equation*}
+   r^\text{grob}(\tilde{v}_h)
+   \colonequals
+   r^\text{fein}(P\tilde{v}_h)
+  \end{equation*}
+  Algebraische Darstellung:
+  \begin{equation*}
+   (r^\text{grob})^T \tilde{x} = (r^\text{fein})^TP\tilde{x}
+   \qquad \implies \qquad
+   r^\text{grob} = P^T r^\text{fein}
+  \end{equation*}
+
+  \item Fehlergleichung auf dem groben Gitter:
+    \begin{equation*}
+     \hat{A}e_c = \hat{r}^k
+     \qquad \text{mit} \qquad
+     \hat{A} = P^T A P
+     \qquad \text{und} \qquad
+     \hat{r}^k = P^T r^k.
+    \end{equation*}
 
-\begin{itemize}
- \item Z.\,B.: Der Jacobi-Vorkonditionierer: $C=\operatorname{diag} (A)^{-1}$
- \item Viele lineare Verfahren werden überhaupt nur betrachtet, um als Vorkonditionierer zu dienen
-       (z.\,B. Mehrgitterverfahren, Gebietszerlegungsverfahren).
 \end{itemize}
-Ein Problem noch: Wie steht es z.\,B.\ mit Gauß--Seidel
-\begin{equation*}
- C = (D-L)^{-1}
- \qquad ?
-\end{equation*}
+
+
+Man erhält die sogenannte \emph{Zweigittermethode}:
+
+\begin{algorithm}[H]
+ \SetAlgoLined
+ % \caption{...}
+ $\hat{A} \longleftarrow P^T A P$
+
+ \For{$k=1,2,3,\dots$}
+ {
+   \textbf{Smooth:} $\nu$ steps of Gauß--Seidel (usually $\nu=3$) to obtain $x^k_*$
+
+   \textbf{Restrict:} compute $\hat{r}^k = P^T r^k = P^T (b - Ax^k_*)$
+
+   \textbf{Coarse correction:} solve $\hat{A}e_c = \hat{r}^k$
+
+   \textbf{Prolong and add:} $x^{k+1} = x^k + Pe_c$
+ }
+\end{algorithm}
+
+\subsection{Mehrgitterverfahren}
+
+Wie lösen wir $\hat{A}e_c = \hat{r}^k$?
+
+Falls $\hat{A}$ klein ist:
 \begin{itemize}
-	\item Funktioniert nicht, denn $W=C^{-1}=D-L$ ist nicht symmetrisch.
+ \item Direkter Löser
 \end{itemize}
-\medskip
-Stattdessen: Symmetrischer Gauß--Seidel
+
+Falls $\hat{A}$ nicht klein ist:
 \begin{itemize}
- \item Abwechselnd
-  \begin{itemize}
-   \item Vorwärtsiteration:
-    \begin{equation*}
-     x^{k+\frac{1}{2}}=x^k+(D-L)^{-1} (b-Ax^k)
-    \end{equation*}
-   %
-   \item Rückwärtsiteration:
-    \begin{equation*}
-     x^{k+1}=x^{k+\frac{1}{2}}+(D-R)^{-1} (b-Ax^{k+\frac{1}{2}})
-    \end{equation*}
-   \end{itemize}
- \item Einsetzen:
-  \begin{equation*}
-   x^{k+1}=x^k+\underbrace{\left[(D-L)^{-1}+(D-R)^{-1}-(D-R)^{-1}A(D-L)^{-1} \right]}_{=C} (b-Ax^k)
-  \end{equation*}
-  Dieses $C$ will man sicher nicht explizit ausrechnen.
- \item Aber: $C^{-1}=W$ ist s.p.d.!
+ \item Mehrgitter!
+ \item Wir müssen $\hat{A}e_c = \hat{r}^k$ nicht exakt lösen!
+ \item Mache rekursiv eine Mehrgitteriteration für $\hat{A}e_c = \hat{r}^k$.
 \end{itemize}
 
+\begin{center}
+ \includegraphics[width=0.5\textwidth]{grid_hierarchy}
+\end{center}
+
+\subsection{Konvergenz}
 
-\section{Glättungseigenschaften}
+\begin{center}
+ \includegraphics[width=0.4\textwidth]{gs_convergence_rates}
+ \includegraphics[width=0.4\textwidth]{mg_convergence_rates}
+
+ Gauß--Seidel \hspace{0.28\textwidth} Multigrid
+\end{center}
 
-\section{Mehrgitterverfahren}
 
 
 \chapter{Teilraumkorrekturverfahren}