Verbessere die Erklärung der Mehrgittermethode

skript-mehrgitter-sander.tex

@@ -18,6 +18,7 @@
 \usepackage{overpic}
 \usepackage{graphicx}
 \usepackage{tikz}
+\usepackage{pgfplots}
 
 \usepackage{amsfonts}
 \usepackage{lmodern}
@@ -953,6 +954,9 @@ Sei $\kappa$ die Kondition der Matrix $A$. Dann gilt nach $k$ Schritten des CG-V
 
 \section{Mehrgitterverfahren: Der anschauliche Zugang}
 
+Das Ziel ist es, einen iterativen Löser zu konstruieren, dessen Konvergenzgeschwindigkeit
+nicht von der Gitterweite $h$ abhängt.
+
 \subsection{Glättung}
 
 Gauß--Seidel-Verfahren und Jacobi-Verfahren sind für typische Finite-Elemente-Matrizen
@@ -962,12 +966,14 @@ extrem langsam.
 
 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
+\smallskip
+
+$\longrightarrow$ Die zwei Verfahren \emph{glätten}.
 
+\bigskip
+
+\emph{Beispiel:} Wir betrachten wieder das Poisson-Problem
 \begin{align*}
  -\Delta u & = 1  \qquad \text{auf $\Omega$}, \\
          u & = 0  \qquad \text{auf $\partial \Omega$,}
@@ -976,8 +982,9 @@ auf dem Einheitsquadrat $\Omega = (0,1)^2$.
 
 \medskip
 
-Hier ist die Finite-Elemente-Lösung auf einem strukturierten Dreiecksgitter
-mit $16 \times 16 \times 2$ Elementen, dargestellt als Höhenplot:
+Das folgende Bild zeigt $u_h$, die Finite-Elemente-Lösung auf einem
+strukturierten Dreiecksgitter mit $16 \times 16 \times 2$ Elementen
+und Lagrange-Elementen erster Ordnung, dargestellt als Höhenplot:
 \begin{center}
  \includegraphics[width=0.6\textwidth,trim=0 50 0 10]{gs_from_zero_iterate50}
 \end{center}
@@ -987,22 +994,6 @@ 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
@@ -1015,17 +1006,17 @@ Hier sind einige der ersten Iterierten:
   \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$ \\
+  $x^0$ & $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$ \\
+  $x^2$ & $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}$
+  $x^{10}$ & $x^{50}$
  \end{tabular}
 \end{center}
 
@@ -1034,7 +1025,7 @@ tatsächlich recht lange.
 
 \bigskip
 
-Jetzt wollen wir zeigen was wir mit \glqq Gauß--Seidel glättet\grqq{} meinen.
+Wir zeigen was wir mit \glqq Gauß--Seidel glättet\grqq{} meinen.
 
 \medskip
 
@@ -1046,17 +1037,17 @@ verrauschten Startwert an.
   \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$ \\
+  $x^0$ & $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$ \\
+  $x^2$ & $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}$
+  $x^{10}$ & $x^{50}$
  \end{tabular}
 \end{center}
 
@@ -1069,8 +1060,26 @@ verrauschten Startwert an.
    verschwunden.
 \end{itemize}
 
+Das kann man sich gut vorstellen: Ein Gauß--Seidel-Schritt ersetzt den Wert $x^k_i$
+an einem Gitterknoten durch
+\begin{equation*}
+ x^{k+1}_i
+  =
+ \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*}
+Dabei tauchen in den Summen nur die Nachbarknoten auf.
+\begin{itemize}
+ \item Zumindest in diesem Modellproblem sind die entsprechenden Koeffizienten
+   alle gleich.
+
+ \item $x^k_i$ wird also durch den Mittelwert seiner Nachbarn ersetzt.
+\end{itemize}
+
+\bigskip
+
+So richtig überzeugend ist diese Argumentation noch nicht:
 
-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?
@@ -1086,65 +1095,84 @@ In Wirklichkeit ist es etwas anders:
  \item Das Gauß--Seidel-Verfahren glättet die \emph{Fehler}!
 \end{itemize}
 
+\missingfigure{Sechs Bilder mit den Fehler für die obigen Iterierten!}
+
 \subsection{Fourier-Analyse}
 \label{sec:jacobi_smoother}
 
-Wir untersuchen jetzt genauer, wie das Jacobi-Verfahren den Fehler glättet.
+Den Glättungseffekt von Gauß--Seidel- und Jacobi-Verfahren kann man besser verstehen,
+wenn man sich wieder die Eigenwerte und Eigenfunktionen von $A$ und $I - CA$ anschaut.
 
-(Das Gauß--Seidel-Verfahren glättet besser, aber es ist auch schwieriger
-zu untersuchen.  Wir lassen es hier aus.)
+\medskip
+
+Wir untersuchen hier das Jacobi-Verfahren.  Für das Gauß--Sidel-Verfahren geht es auch,
+ist aber etwas komplizierter.  Interessierte finden die Rechnung zum Beispiel bei
+\citet{trottenberg}
+\todo[inline]{Kapitelnummer!}
 
 Zur Erinnerung: Das gedämpfte Jacobi-Verfahren ist
 \begin{equation*}
  x^{k+1}
  =
- (I - \omega D^{-1}A)x^k + D^{-1}b.
+ (I - \eta D^{-1}A)x^k + D^{-1}b.
 \end{equation*}
 
 Wir betrachten wieder das Modellproblem.  Dort ist
 \begin{equation*}
- A_ii = \frac{4}{h^2}
+ D_i \colonequals A_{ii} = \frac{4}{h^2}
  \qquad
  \text{für alle $i=1,\dots,n$},
 \end{equation*}
 und somit
 \begin{equation*}
- x^{k+1} = (I- \omega \frac{h^2}{4} A) x^k + \frac{h^2}{4}b.
+ x^{k+1} = (I- \eta \frac{h^2}{4} A) x^k + \frac{h^2}{4}b.
 \end{equation*}
 
 Um das Glättungsverhalten des Jacobi-Verfahrens zu verstehen berechnen wir
 die Eigenfunktionen der Iterationsmatrix $I - \omega \frac{h^2}{4}A$.
 
+\medskip
+
 Die Eigenfunktionen von $I - \omega \frac{h^2}{4}A$ sind offensichtlich
 die selben wie die von $A$, und zwar
 \begin{equation*}
- \varphi_h^{i,j}
+ \varphi_h^{\nu,\mu}
  =
- \sin i\pi x \cdot \sin j \pi y.
+ \sin \nu\pi x \cdot \sin \mu \pi y.
 \end{equation*}
 Die Eigenwerte von $A$ sind
-\todo[inline]{Fehlen noch}
-Die Eigenwerte von $I - \omega \frac{h^2}{4}A$ sind
 \begin{equation*}
- \chi_h^{i,j}
+ \lambda{\nu\mu}
  =
- 1 - \frac{\omega}{2} ( 2 - \cos i \pi h - \cos j \pi h)
+ \frac{4}{h^2} \Big[ \sin^2\big(\frac{1}{2} \nu \pi h \big) + \sin^2\big(\frac{1}{2} \mu \pi h \big)\Big]
+ \qquad
+ \forall 1 \le \nu,\mu < n.
 \end{equation*}
 
+Die Eigenwerte von $I - \eta \frac{h^2}{4}A$ sind
+\begin{equation*}
+ \chi_h^{\nu,\mu}
+ =
+ 1 - \eta \frac{h^2}{4} \lambda^{\nu,\mu}
+ =
+ 1 - \frac{\eta}{2} ( 2 - \cos \nu \pi h - \cos \mu \pi h)
+\end{equation*}
+(wegen der Rechenregel $2\sin^2 \alpha = 1 - \cos 2\alpha$).
+
 \bigskip
 
 Die Konvergenzgeschwindigkeit des Verfahrens wird bekanntlich durch den
-Spektralradius von $I-CA$ bestimmt, und der ist hier für alle $0 < \omega \le 1$
+Spektralradius von $I-CA$ bestimmt, und der ist hier für alle $0 < \eta \le 1$
 \begin{equation*}
  \rho (I-CA)
  =
  \abs{\chi_h^{1,1}}
  =
- \abs{1 - \omega(1-\cos \pi h)}
+ \abs{1 - \eta(1-\cos \pi h)}
  =
- 1 - O(\omega h^2).
+ 1 - O(\eta h^2).
 \end{equation*}
-Falls $\omega \le 0$ oder $\omega >1$ gilt
+Falls $\eta \le 0$ oder $\eta >1$ gilt
 \begin{equation*}
  \rho(I-CA) > 1
 \end{equation*}
@@ -1154,41 +1182,133 @@ Falls $\omega \le 0$ oder $\omega >1$ gilt
 
 Jetzt schreiben wir die Fehlerfortpflanzung
 \begin{equation*}
- e^{k+1} = x^{k+1} - x^*
+ e^{k+1} \colonequals x^* - x^{k+1}
  =
- (I - \omega \frac{h^2}{4}A)(x^k - x^*)
+ \Big(I - \eta \frac{h^2}{4}A\Big)(x^* - x^k)
 \end{equation*}
-in der Eigenvektorbasis
+in der Eigenvektorbasis der Iterationsmatrix
 \begin{align*}
  e^k
  & =
- \sum_{i,j = 1}^n \alpha_{i,j}^k \cdot \varphi_h^{i,j} \\
+ \sum_{\nu,\mu = 1}^n \alpha_{\nu,\mu}^k \cdot \varphi_h^{\nu,\mu} \\
  %
  e^{k+1}
  & =
- (I - \omega \frac{h^2}{4}A) e^k
+ \Big(I - \eta \frac{h^2}{4}A\Big) e^k
  =
- \sum_{i,j = 1}^n \chi_h^{i,j} \alpha_{i,j}^k \cdot \varphi_h^{i,j}.
+ \sum_{\nu,\mu = 1}^n \chi_h^{\nu,\mu} \alpha_{\nu,\mu}^k \cdot \varphi_h^{\nu,\mu}.
 \end{align*}
 
-Die Eigenfunktionen sind \glqq hochfrequent\grqq{} für große $i$, $j$,
-und \glqq niederfrequent\grqq{} für kleine $i$, $j$.
-
-\medskip
-
-Gleichzeitig sind die Eigenwerte für große $i$, $j$ klein, und für kleine $i$, $j$, groß.
-\missingfigure{Eigenwerte als Funktion von $i$ und $j$.}
-
-Deshalb verschwinden die hochfrequenten Anteile des Fehlers schnell,
-obwohl der Fehler insgesamt nur sehr langsam verschwindet.
-
+Die Eigenfunktionen sind \glqq hochfrequent\grqq{} für große $\nu$, $\mu$,
+und \glqq niederfrequent\grqq{} für kleine $\nu$, $\mu$.
 
+\begin{center}
+ \begin{tabular}{ccc}
+ \begin{tikzpicture}[scale=0.5]
+  \begin{axis}[zmin=-1,zmax=1]
+   \addplot3 [surf,
+              samples=16,
+              domain=0:1,y domain=0:1
+             ] {sin(deg(pi * x)) * sin(deg(pi*y))};
+  \end{axis}
+ \end{tikzpicture}
+ &
+ \begin{tikzpicture}[scale=0.5]
+  \begin{axis}
+   \addplot3 [surf,
+              samples=16,
+              domain=0:1,y domain=0:1
+             ] {sin(deg(pi * x)) * sin(deg(2*pi*y))};
+  \end{axis}
+ \end{tikzpicture}
+ &
+ \begin{tikzpicture}[scale=0.5]
+  \begin{axis}
+   \addplot3 [surf,
+              samples=16,
+              domain=0:1,y domain=0:1
+             ] {sin(deg(pi * x)) * sin(deg(3*pi*y))};
+  \end{axis}
+ \end{tikzpicture}
+ \\
+ \begin{tikzpicture}[scale=0.5]
+  \begin{axis}
+   \addplot3 [surf,
+              samples=16,
+              domain=0:1,y domain=0:1
+             ] {sin(deg(2*pi * x)) * sin(deg(pi*y))};
+  \end{axis}
+ \end{tikzpicture}
+ &
+ \begin{tikzpicture}[scale=0.5]
+  \begin{axis}
+   \addplot3 [surf,
+              samples=16,
+              domain=0:1,y domain=0:1
+             ] {sin(deg(2*pi * x)) * sin(deg(2*pi*y))};
+  \end{axis}
+ \end{tikzpicture}
+ &
+ \begin{tikzpicture}[scale=0.5]
+  \begin{axis}
+   \addplot3 [surf,
+              samples=16,
+              domain=0:1,y domain=0:1
+             ] {sin(deg(2*pi * x)) * sin(deg(3*pi*y))};
+  \end{axis}
+ \end{tikzpicture}
+ \\
+ \begin{tikzpicture}[scale=0.5]
+  \begin{axis}
+   \addplot3 [surf,
+              samples=16,
+              domain=0:1,y domain=0:1
+             ] {sin(deg(3*pi * x)) * sin(deg(pi*y))};
+  \end{axis}
+ \end{tikzpicture}
+ &
+ \begin{tikzpicture}[scale=0.5]
+  \begin{axis}
+   \addplot3 [surf,
+              samples=16,
+              domain=0:1,y domain=0:1
+             ] {sin(deg(3*pi * x)) * sin(deg(2*pi*y))};
+  \end{axis}
+ \end{tikzpicture}
+ &
+ \begin{tikzpicture}[scale=0.5]
+  \begin{axis}
+   \addplot3 [surf,
+              samples=16,
+              domain=0:1,y domain=0:1
+             ] {sin(deg(3*pi * x)) * sin(deg(3*pi*y))};
+  \end{axis}
+ \end{tikzpicture}
+ \end{tabular}
+\end{center}
 
+\medskip
 
+Gleichzeitig sind die Eigenwerte für große $\nu$, $\mu$ klein, und für kleine $\nu$, $\mu$, groß.
 
+\todo[inline]{Noch mal nachprüfen -- das stimmt gar nicht...}
 
+Der folgende Plot zeigt sie für $\eta = 1$:
+\begin{center}
+ \begin{tikzpicture}
+  \begin{axis}[colorbar]
+   \addplot3 [surf,
+%              faceted color=blue,
+              samples=15,
+              domain=1:15,y domain=1:15
+             ] {1 - 0.5*(2 - cos (deg(x*pi/16)) - cos (deg(y*pi/16)))};
+  \end{axis}
+ \end{tikzpicture}
+\end{center}
 
 
+Deshalb verschwinden die hochfrequenten Anteile des Fehlers schnell,
+obwohl der Fehler insgesamt nur sehr langsam verschwindet.
 
 
 \subsection{Zweigitterverfahren}
@@ -1196,17 +1316,22 @@ obwohl der Fehler insgesamt nur sehr langsam verschwindet.
 Die zweite wichtige Einsicht ist:
 
 \begin{itemize}
- \item Der Fehler $e^k \colonequals x^* - x^k$ löst eine lineare Gleichung
+ \item Angenommen wir haben als aktuelle Iterierte $x^k$.  Statt der Lösung $x^*$
+  könnten wir auch versuchen, den aktuellen Fehler $e^k \colonequals x^* - x^k$
+  auszurechnen.
+
+ \item Der Fehler $e^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.
+  \item Da $e^k$ glatt ist kann man eine gute Approximation erwarten, wenn man
+   die Fehlergleichung $Ae = r^k$ auf einem gröberen Gitter löst.
 \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.
+von $V_h^\text{fein} = V_h$ ist.
 
    \begin{center}
     \begin{tikzpicture}[scale=0.35]
@@ -1230,10 +1355,14 @@ von $V_h^\text{fein} = V_H$ ist.
    \end{center}
 
 \begin{itemize}
-  \item Definiere Prolongationsoperator als die kanonische Injektion
+ \item Wir definieren eine Abbildung aus dem groben in den feinen Funktionenraum
     \begin{equation*}
-     P : V_h^\text{grob} \to V_h^\text{fein}
+     P : V_h^\text{grob} \to V_h^\text{fein}.
     \end{equation*}
+  Man nennt $P$ den \emph{Prolongationsoperator}.
+
+ \item Da in diesem Beispiel die Räume geschachtelt sind nehmen wir einfach
+  die kanonische Injektion:
 
     \begin{center}
      \includegraphics[width=0.3\textwidth]{coarse_function}
@@ -1242,37 +1371,43 @@ von $V_h^\text{fein} = V_H$ ist.
     \end{center}
 \end{itemize}
 
-Mit $P$ kann man Funktionen vom groben auf das feine Gitter transferieren,
+Wenn $P$ wie hier die kanonische Injektion ist kann man Funktionen vom groben auf das feine Gitter transferieren,
     ohne einen Fehler zu machen.
 
-\medskip
+\subsubsection{Konstruktion des Grobgitterproblems}
 
+Um die Fehlergleichung auf dem groben Gitter lösen zu können brauchen wir
+eine algebraische Darstellung.
 \begin{itemize}
+ \item Genauer: Wir brauchen eine Matrix $A^H$ und einen Vektor $r^H$,
+  der der Bilinearform $a(\cdot,\cdot)$ und dem Residuum $\ell(\cdot) - a(u_h^k,\cdot)$
+  auf dem groben Gitter entspricht.
+
  \item Natürlich können wir keine Funktion vom feinen Gitter auf das grobe Gitter
   übertragen, ohne einen Approximationsfehler zu machen.
 
- \item \emph{Aber:} Wir können Linearformen $\ell(\cdot) : V_h^\text{fein} \to \R$
+ \item \emph{Aber:} Wir können Linearformen $\mathcal{L}(\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)
+   \mathcal{L}^\text{grob}(v_h) \colonequals \mathcal{L}^\text{fein}(Pv_h)
    \qquad
    \forall v_h \in V_h^\text{grob}.
   \end{equation*}
 
- \item Das Residuum $r^k$ ist die algebraische Darstellung so einer Linearform:
+ \item Das Residuum $r^\text{fein} \colonequals r^k$ ist die algebraische Darstellung so einer Linearform:
 
-  Funktionenraum:
+  Im Funktionenraum:
   \begin{equation*}
-   r^k(v_h) \colonequals \langle b,v_h \rangle - a(u_h^k,v_h)
+   r^k(v_h) \colonequals \ell(v_h) - a(u_h^k,v_h)
    \qquad
    \forall v_h \in V_h^\text{fein}
   \end{equation*}
 
   Algebraisch :
   \begin{equation*}
-   (r^k)^T v \colonequals b^T v - (x^k)^T A v
+   (r^k)^T y \colonequals b^T y - (x^k)^T A y
    \qquad
-   \forall v \in \R^n
+   \forall y \in \R^{n_\text{fein}}.
   \end{equation*}
 
  \item Die Einschränkung auf den groben Funktionenraum ist
@@ -1280,12 +1415,17 @@ Mit $P$ kann man Funktionen vom groben auf das feine Gitter transferieren,
    r^\text{grob}(\tilde{v}_h)
    \colonequals
    r^\text{fein}(P\tilde{v}_h)
+   \qquad
+   \forall \tilde{v}_h \in V_h^\text{grob}.
   \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}
+   (r^\text{grob})^T \tilde{y} = (r^\text{fein})^TP\tilde{y}
+   \qquad \forall \tilde{y} \in \R^{n_\text{grob}},
+  \end{equation*}
+  und deshalb
+  \begin{equation*}
+   r^\text{grob} = P^T r^\text{fein}.
   \end{equation*}
 
   \item Fehlergleichung auf dem groben Gitter:
@@ -1299,6 +1439,18 @@ Mit $P$ kann man Funktionen vom groben auf das feine Gitter transferieren,
 
 \end{itemize}
 
+Die Abbildung $A \mapsto P^T A P$ heißt \emph{Galerkin-Restriktion}.
+
+\begin{exercise}
+ Zeigen Sie dass $P^TAP$ positiv definit ist wenn $A$ positiv definit ist.
+\end{exercise}
+
+\begin{exercise}€
+ Programmieren Sie das Produkt $P^TAP$, wenn $P$ und $A$ im CRS-Format
+ gegeben sind.
+\end{exercise}
+
+
 
 Man erhält die sogenannte \emph{Zweigittermethode}: