skript-mehrgitter-sander.tex 160 KB
Newer Older
1
% The following line selects the page layout for printing
2
\documentclass{scrbook}
3
4
5

% The following lines select a page layout without vertical margins,
% and a very large right margin.  This is much better for online teaching.
6
7
8
9
% \documentclass[paper=310mm:297mm]{scrbook}
% \KOMAoptions{twoside=false}
%
% \usepackage[width=12cm, left=3cm, top=0.2cm, bottom=0.5cm]{geometry}
10
11
12
13

\usepackage[utf8]{inputenc}

\usepackage{longtable,tabularx,tabulary,booktabs}
14
\usepackage{aligned-overset}
15
\usepackage{array}
16
\usepackage[linesnumbered,inoutnumbered]{algorithm2e}
17
18
\usepackage{caption}
\usepackage{url}
Sander, Oliver's avatar
Sander, Oliver committed
19
\usepackage{overpic}
20
21
\usepackage{graphicx}
\usepackage{tikz}
22
\usepackage{pgfplots}
Sander, Oliver's avatar
Sander, Oliver committed
23
\pgfplotsset{compat=1.16}
24
25
26
27

\usepackage{amsfonts}
\usepackage{lmodern}
\usepackage[ngerman]{babel}
28
\usepackage{csquotes}
29
30
31
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
32
\usepackage{bm}
33
\usepackage{mathtools}
34
\usepackage{mathabx}   % for \vvvert
35
36
37
\usepackage{colonequals}
\usepackage{enumitem}

38
\usepackage[maxbibnames=99,  % Never write 'et al.' in a bibliography
39
40
41
42
43
44
45
            giveninits=true, % Abbreviate first names
            sortcites=true,  % Sort citations when there are more than one
            natbib=true,     % Allow commands like \citet and \citep
            backend=biber]{biblatex}
\addbibresource{skript-mehrgitter-sander.bib}

\usepackage{microtype}
46
\usepackage[disable]{todonotes}
47
48
49
50
51
52
53
54
55
56
\usepackage{esdiff}
\usepackage{hyperref}

\allowdisplaybreaks[0]
\setlength{\parindent}{0pt}  % No indent of first line of paragraph

\newcommand{\R}{\mathbb R}
\newcommand{\N}{\mathbb N}
\newcommand{\curl}{\operatorname{curl}}
\renewcommand{\div}{\operatorname{div}}
57
\newcommand{\grad}{\operatorname{grad}}
58
59
\newcommand{\tr}{\operatorname{tr}}
\DeclareMathOperator*{\argmin}{arg\,min}
Sander, Oliver's avatar
Sander, Oliver committed
60
\newcommand{\Gammatight}[1]{{\Gamma\hspace{-0.8mm}_{#1}}}  % A \Gamma with a subscript, and extra kerning
61

62
63
64
\newcommand{\Radd}{R^\textup{add}}
\newcommand{\Rmult}{R^\textup{mult}}

65
% Bold letters
66
\newcommand{\bj}{\bm{j}}
67
\newcommand{\bn}{\mathbf{n}}
68
\newcommand*\bu{\bm{u}}
69
70
\newcommand{\bv}{\mathbf{v}}

71
\newcommand*\bA{\bm{A}}
72
73
74
75
\newcommand*\bB{\bm{B}}
\newcommand*\bD{\bm{D}}
\newcommand*\bE{\bm{E}}
\newcommand*\bF{\bm{F}}
76
77
78
79
80
81
\newcommand*\bH{\bm{H}}
\newcommand*\bI{\bm{I}}
\newcommand*\bL{\bm{L}}
\newcommand*\bP{\bm{P}}
\newcommand*\bQ{\bm{Q}}
\newcommand*\bV{\bm{V}}
82
83
84

\newcommand{\btau}{\bm{\tau}}

85
86
87
88
89
90
91
92
93
94
95
96
\newcommand*\bHd{\bH(\div)}
\newcommand*\til{\widetilde}
\newcommand{\stackrela}[2]{\overset{#1}&{#2}}

\newcommand{\X}{\mathcal{X}}				 % Triangulierung
\newcommand{\V}{\mathcal{V}}				 % Triangulierung
\newcommand{\E}{\mathcal{E}}				 % Triangulierung
\newcommand{\F}{\mathcal{F}}				 % Triangulierung
\newcommand{\T}{\mathcal{T}}				 % Triangulierung

\DeclarePairedDelimiter{\set}{\lbrace}{\rbrace}
\DeclarePairedDelimiter{\norma}{\Vert}{\Vert_{A}}
97
98
\DeclarePairedDelimiter{\abs}{\lvert}{\rvert}
\DeclarePairedDelimiter{\norm}{\lVert}{\rVert}
99
\DeclarePairedDelimiter{\triplenorm}{\vvvert}{\vvvert}
100
101
102

\theoremstyle{plain} %Text ist Kursiv
\newtheorem{theorem}{Satz}[chapter]
103
\newtheorem{hilfssatz}[theorem]{Hilfssatz}
104
105
106
107
108
109
110
111
112
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Korollar}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{problem}[theorem]{Problem}

\theoremstyle{definition} %Text ist \"upright"
\newtheorem{remark}[theorem]{Bemerkung}
\newtheorem{example}[theorem]{Beispiel}
113
\newtheorem{exercise}[theorem]{Übung}
114
115
116
117
118
119
120
121
122
123
124
125
126
127

\swapnumbers

% Suche auch in diesem Verzeichnis nach Bildern
\graphicspath{{gfx/}}

\begin{document}

\begin{titlepage}

\begin{center}

\vspace*{0.3\textheight}

128
{\Huge Mehrgitterverfahren}
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
\vspace*{0.05\textheight}

{\Large Oliver Sander}
\vspace*{0.1\textheight}

\today
\vspace{0.4\textheight}

\end{center}

\end{titlepage}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Copyright notice
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\thispagestyle{empty}

\vspace*{\fill}
\noindent
\includegraphics[width=0.2\textwidth]{license-logo-cc-by-sa}
Oliver Sander, 2021

\bigskip

\noindent
Copyright 2021 by Oliver Sander.
This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License.
To view a copy of this license, visit \url{http://creativecommons.org/licenses/by-sa/4.0/}.

\pagebreak

\tableofcontents
\thispagestyle{empty}
\setcounter{page}{0}

Sander, Oliver's avatar
Sander, Oliver committed
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
\chapter{Finite Elemente Methoden}

In diesem Kapitel geben wir eine kurze Einführung in die Methode der Finiten Elemente.
Damit legen wir die Grundlagen, und motivieren die weiteren Auführungen.

\section{Schwache Formulierung}

\begin{figure}
 \begin{center}
  \begin{overpic}[width=0.3\textwidth]{gfx/gebiet}
   \put(50,82){$\Gamma_D$}
   \put(25,2){$\Gamma_N$}
   \put(30,40){$\Omega$}
  \end{overpic}
 \end{center}
\caption{Einfaches Gebiet $\Omega$ mit Dirichletrand $\Gamma_D$ und Neumannrand
$\Gamma_N$.}
\label{fig:gebiet}
\end{figure}

Sei $\Omega \subset \R^d$ ein Gebiet und $f : \Omega \to \R$.
Der Rand $\partial \Omega$ von $\Omega$ sei zerlegt in zwei disjunkte Teile $\Gammatight{D}$ und $\Gammatight{N}$,
und wir bezeichnen die äußere Einheitsnormale mit $\nu$.
187
Wir betrachten als Beispiel die \emph{Poisson-Gleichung}
Sander, Oliver's avatar
Sander, Oliver committed
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
\begin{equation}
\label{eq:starke_formulierung}
 - \Delta u \colonequals - \sum_{i=1}^d \frac{\partial^2 u}{\partial x_i^2} = f
\qquad \text{auf $\Omega$},
\end{equation}
mit den Randbedingungen
\begin{align}
\nonumber
 u & = 0 \qquad \text{auf $\Gammatight{D}$}, \\
 \label{eq:neumann_randbedingungen}
\frac{\partial u}{\partial \nu} & = g \qquad \text{auf $\Gammatight{N}$}.
\end{align}
Wir bezeichnen mit $H^1(\Omega)$ den Raum aller skalarer $L^2$-Funktionen auf
$\Omega$ mit schwacher Ableitung in $L^2(\Omega)$.
Als Raum von Testfunktionen führen wir
\begin{equation*}
 H^1_{\Gammatight{D}} (\Omega)
=
\{ v \in H^1(\Omega) \; | \; v|_{\Gammatight{D}} = 0 \}
\end{equation*}
ein, wobei $v|_{\Gammatight{D}}$ im Sinne des Spursatzes zu verstehen ist.
Zur Abkürzung setzen wir $H \colonequals H^1_{\Gammatight{D}} (\Omega)$.

Multiplikation von \eqref{eq:starke_formulierung} mit einer Testfunktion
$v \in H^1_{\Gammatight{D}} (\Omega)$ und Integration über $\Omega$ ergibt
\begin{equation*}
 - \int_\Omega \Delta u \cdot v \, dx = \int_\Omega f v \, dx.
\end{equation*}
Wir wenden darauf die Greensche Formel an, und erhalten
\begin{equation*}
 \int_\Omega \nabla u \nabla v \, dx
=
\int_{\partial \Omega} \frac{\partial u}{\partial \nu} v \, ds
+
\int_\Omega f v \, dx.
\end{equation*}
Nach Voraussetzung besteht $\partial \Omega$ aus den zwei Teilen $\Gammatight{D}$ und $\Gammatight{N}$.
Die Testfunktion $v$ wurde so konstruiert, dass $v|_{\Gammatight{D}} = 0$.  Mithin bleibt
nur ein Integral über $\Gammatight{N}$ übrig.  Dort können wir die Neumannrandbedingungen~\eqref{eq:neumann_randbedingungen}
einsetzen und erhalten
\begin{equation*}
 \int_\Omega \nabla u \nabla v \, dx
=
\int_{\Gamma_N} g v \, ds
+
\int_\Omega f v \, dx.
\end{equation*}
Wir schreiben dies als
\begin{equation}
\label{eq:schwache_formulierung}
238
 a(u,v) = \ell(v)
Sander, Oliver's avatar
Sander, Oliver committed
239
240
241
\end{equation}
mit
\begin{equation*}
242
243
244
 a(v,w) \colonequals \int_\Omega \nabla v \nabla w \, dx
\qquad \text{und} \qquad
 \ell(v) \colonequals \int_{\Gamma_N} g v \, ds
Sander, Oliver's avatar
Sander, Oliver committed
245
246
247
+
\int_\Omega f v \, dx.
\end{equation*}
248
249
Man nennt \eqref{eq:schwache_formulierung} die \emph{schwache Formulierung}
(oder \emph{Variationsformulierung})
Sander, Oliver's avatar
Sander, Oliver committed
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
von \eqref{eq:starke_formulierung}.

$a(\cdot, \cdot)$ ist eine symmetrische Bilinearform.  Hat $\Gammatight{D}$
positives $d-1$-dimensionales Maß, so ist $a(\cdot,\cdot)$
$H^1_{\Gamma_D}$-elliptisch (siehe z.B.\ \cite{braess:2013}), d.h.\ es existieren $\gamma, \Gamma > 0$
so daß
\begin{equation*}
 \gamma \norm{v}^2_H \le a(v,v)
\qquad
\abs{a(v,w)} \le \Gamma \norm{v}_H \norm{w}_H
\qquad
\text{für alle $v,w \in H^1_{\Gamma_D}(\Omega)$}.
\end{equation*}


\begin{theorem}[Lax-Milgram-Lemma \cite{braess:2013}]
Sei $a(\cdot,\cdot)$ $H$-elliptisch (nicht notwendig symmetrisch!).
Dann hat die Variationsgleichung
\begin{equation*}
 u \in H
\quad : \quad
271
a(u,v) = \ell(v)
Sander, Oliver's avatar
Sander, Oliver committed
272
273
\qquad v \in H
\end{equation*}
274
for jedes $\ell \in H'$ eine eindeutig bestimmte Lösung.
Sander, Oliver's avatar
Sander, Oliver committed
275
276
277
278
279
280
281
282
283
\end{theorem}




\section{Diskretisierung (Galerkin-Verfahren)}

Sei $V_h$ ein endlichdimensionaler Teilraum von $H^1_{\Gammatight{D}}$ (mit Dimension $n$).
Der Index $h$
284
bezeichne einen Ordnungsparameter.  Die \emph{diskrete Variationsformulierung}
Sander, Oliver's avatar
Sander, Oliver committed
285
286
287
lautet
\begin{equation}
\label{eq:diskrete_formulierung}
288
 u_h \in V_h \quad : \quad a(u_h,v_h) = \ell(v_h)
Sander, Oliver's avatar
Sander, Oliver committed
289
290
291
292
293
294
295
296
297
\qquad \text{für alle $v_h \in V_h$}.
\end{equation}
Anwendung des Lax--Milgram-Lemmas ergibt wieder Existenz und Eindeutigkeit
von Lösungen.

Sei $\{\varphi_i\}, 0 \le i < n$, eine Basis von $V_h$. Dann ist
\eqref{eq:diskrete_formulierung} äquivalent zu
\begin{equation*}
%\label{eq:diskrete_formulierung}
298
 u_h \in V_h \quad : \quad a(u_h, \varphi_i) = \ell(\varphi_i)
Sander, Oliver's avatar
Sander, Oliver committed
299
300
301
\qquad \text{für alle $0 \le i < n$}.
\end{equation*}
Dies wiederum entspricht dem linearen Gleichungssystem
302
303
\begin{equation}
\label{eq:fe_lineares_gleichungssystem}
304
 Ax = b,
305
\end{equation}
Sander, Oliver's avatar
Sander, Oliver committed
306
307
308
wobei
\begin{alignat*}{2}
A &\in \R^{n \times n},
309
\qquad &A_{ij} & \colonequals a(\varphi_i, \varphi_j)
Sander, Oliver's avatar
Sander, Oliver committed
310
311
312
313
314
=
\int_\Omega \nabla \varphi_i \nabla \varphi_j \, dx \\
%
b &\in \R^n,
\qquad &
315
b_i & \colonequals \ell(\varphi_i)
Sander, Oliver's avatar
Sander, Oliver committed
316
317
318
319
320
=
\int_{\Gamma_N} g \varphi_i \, ds
+
\int_\Omega f \varphi_i \, dx
\end{alignat*}
321
und $x \in \R^n$ die Koeffizienten von $u_h$ bzgl.\ der Basis $\{\varphi_i\}$ sind.
Sander, Oliver's avatar
Sander, Oliver committed
322

323
324
325
\bigskip

Die Matrix $A$ nennt man \emph{Steifigkeitsmatrix}.  Dieser Ausdruck kommt aus den
Sander, Oliver's avatar
Sander, Oliver committed
326
327
328
Anfangstagen der Finiten Elemente in der Ingenieurswelt, wo insbesondere mechanische
Probleme gelöst wurden.  Die Matrix $A$ beschreibt dabei die Steifigkeit des simulierten
Objekts, und der Vektor $b$ beschreibt die von außen wirkenden Kräfte.  Den Vektor $b$
329
nennt man deshalb manchmal auch \emph{Lastvektor}.
Sander, Oliver's avatar
Sander, Oliver committed
330
331
332
333
334
335
336
337


\section{Finite Elemente}
\label{sec:finite_elemente}

Finite Elemente sind eine bestimmte Möglichkeit, den Teilraum
$V_h$ zu wählen.

338
339
340
341
\medskip

Sei $\Omega$ von jetzt ab \emph{polygonal berandet}.  Wir füllen
$\Omega$ mit einer Triangulierung / einem \emph{Gitter}.  Ein Beispielgitter für
Sander, Oliver's avatar
Sander, Oliver committed
342
343
344
345
346
347
348
349
350
351
352
353
354
355
ein zweidimensionales Gebiet sieht man in Abbildung~\ref{fig:2d-Beispielgitter}.

\begin{figure}
\begin{center}
\includegraphics[width=0.5\textwidth]{dreiecksgitter}
\end{center}
\caption{Ein einfaches zweidimensionales Gitter}
\label{fig:2d-Beispielgitter}
\end{figure}

\begin{definition}
\label{def:triangulierung}
 Es sei $\Omega \subset \R^2$ polygonal berandet, und $\Gammatight{D}$ auf dem
Rand aufgelöst.  Eine Menge $\mathcal{T}$
356
von abgeschlossenen Dreiecken $T$ heißt \emph{Triangulierung} von $\Omega$, falls gilt
Sander, Oliver's avatar
Sander, Oliver committed
357
358
359
\begin{enumerate}
 \item Die Menge der Dreiecke überdeckt den Abschluss des Gebiets
\begin{equation*}
360
 \overline{\Omega} = \bigcup_{T \in \mathcal{T}} T.
Sander, Oliver's avatar
Sander, Oliver committed
361
362
363
364
365
366
367
368
\end{equation*}
%
\item Der Schnitt zweier Dreiecke aus $\mathcal T$ ist entweder eine
gemeinsame Kante, ein gemeinsamer Eckpunkt, oder leer.
\end{enumerate}
\end{definition}
Der Diskretisierungsparameter
\begin{equation*}
369
 h \colonequals \max_{T \in \mathcal{T}} \operatorname{diam} T,
Sander, Oliver's avatar
Sander, Oliver committed
370
371
\end{equation*}
der als Index in $V_h$ auftaucht, beschreibt traditionell die Feinheit des Gitters.
372
Oft nennt man $h$ auch \emph{Gitterweite}.
Sander, Oliver's avatar
Sander, Oliver committed
373
374
375
376

\begin{definition}
Der Raum
\begin{equation*}
377
 V_h^{p}
Sander, Oliver's avatar
Sander, Oliver committed
378
379
380
\colonequals
\{ v \in C(\overline{\Omega})
\; | \;
381
\text{$v$ ist Polynom von Grad $\le p$ auf $T$ für alle $T \in \mathcal{T}$, und $v|_{\Gamma_D} = 0$} \}
Sander, Oliver's avatar
Sander, Oliver committed
382
\end{equation*}
383
heißt Raum der Lagrange-Elemente $p$-Ordnung
Sander, Oliver's avatar
Sander, Oliver committed
384
385
386
387
388
389
390
391
bezüglich der Triangulierung $\mathcal{T}$.
\end{definition}

\begin{figure}
\begin{center}
\includegraphics[width=0.4\textwidth]{generalP1function}
\includegraphics[width=0.4\textwidth]{nodalP1basis}
\end{center}
392
\caption{Links: eine Finite-Elemente Funktion.  Rechts: ein Element der Knotenbasis erster Ordnung}
Sander, Oliver's avatar
Sander, Oliver committed
393
394
395
\label{fig:beispiel-fe-funktion}
\end{figure}

396
397
398
Als Basis von $V_h^{p}$ wählt man gewöhnlich die \emph{Knotenbasis}.  Sei $\mathcal L$
die Menge der Lagrange-Knoten des Gitters.  Dann ist die Knotenbasis die Menge
der Funktionen $\varphi_i \in V_h^{p}$ mit
Sander, Oliver's avatar
Sander, Oliver committed
399
400
401
\begin{equation*}
 \varphi_i (v_j) = \delta_{ij}
\qquad
402
\text{für alle $v_j \in \mathcal{L}$}.
Sander, Oliver's avatar
Sander, Oliver committed
403
\end{equation*}
404
Es ist also insbesondere $\dim V_h^{p} = \abs{\mathcal L \setminus \mathcal{L}_{\Gamma_D}}$
Sander, Oliver's avatar
Sander, Oliver committed
405
406
(Abbildung~\ref{fig:beispiel-fe-funktion}).

407
408
409
410
\bigskip

Der Diskretisierungsfehler hängt hauptsächlich von der Gitterauflösung $h$,
der Ordnung $p$ der Finite-Elemente-Funktionen, und der Glattheit von $u$ ab.
Sander, Oliver's avatar
Sander, Oliver committed
411
412

\begin{theorem}
413
414
415
 Es sei $h$ klein genug, $u \in H^1_0(\Omega) \cap H^{p+1}(\Omega)$ die Lösung
 der schwachen Formulierung, und $u_h$ die entsprechende Finite-Elemente-Lösung in $V_h^p$.
 Dann gilt die a priori Abschätzung
Sander, Oliver's avatar
Sander, Oliver committed
416
\begin{equation*}
417
 \norm{u - u_h}_1 \le C h^p \abs{u}_{p+1}
Sander, Oliver's avatar
Sander, Oliver committed
418
419
420
421
422
423
424
425
\end{equation*}
mit einer Konstanten $C$.
\end{theorem}

\begin{proof}
 \citet{braess:2013}
\end{proof}

426
\chapter{Mehrgitterverfahren}
Sander, Oliver's avatar
Sander, Oliver committed
427

428
\section{Iterative Lösungsverfahren für dünnbesetzte lineare Gleichungssysteme}
429
\label{sec:iterative_loeser}
430

431
432
Wir beschäftigen uns jetzt mit der Frage, wie die linearen
Gleichungssysteme~\eqref{eq:fe_lineares_gleichungssystem} gelöst werden können.
433

434
\subsection{Eigenschaften von Finite-Elemente-Matrizen}
435
436
437

Finite-Elemente-Matrizen haben bestimmte Eigenschaften, die das Lösen schwierig machen.

438
\subsubsection{Größe}
439
Finite-Elemente-Matrizen können sehr groß werden:
440
\begin{itemize}
441
442
443
444
445
446
447
 \item Für jeden Lagrangeknoten enthält~\eqref{eq:fe_lineares_gleichungssystem} eine Gleichung.

 \item Für $d$-dimensionale Gebiete und Lagrange-Elemente erster Ordnung
       hat man etwa $n \approx \operatorname{Vol}( \Omega ) \cdot h^{-d}$ Knoten.

 \item Je feiner das Gitter, desto präziser die Lösung, desto größer aber auch die Matrix.

448
 \item Eine Zahl in doppelter Genauigkeit braucht $8$~Byte.  In ein Gigabyte Hauptspeicher
449
450
451
452
453
454
455
456
457
   passen also ca.\ $1{,}25\cdot 10^{8}$ Zahlen.
\end{itemize}

\begin{exercise}
 Wie hängt die Problemgröße von der Ordnung der Ansatzfunktionen ab?
\end{exercise}

\subsubsection{Kondition}

458
459
Die Matrizen zu elliptischen Problemen sind positiv definit, aber sie
sind schlecht konditioniert.
460

461
\medskip
462

463
464
465
466
467
468
469
470
471
472
473
474
475
\emph{Beispiel:} Die Poisson-Gleichung auf $\Omega = [0,1]^2$, Finite-Elemente
mit Lagrange-Elementen auf einem strukturierten Dreiecksgitter mit Kantenlänge $h = \frac{1}{n}$.
\begin{itemize}
 \item Die Eigenwerte der Steifigkeitsmatrix $A$ sind~\cite[Kapitel~12.3.3]{dahmen_reusken:2008}
  \begin{align*}
   \lambda_{\nu \mu}
   & =
   \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
   1 \le \nu,\mu < n.
  \end{align*}

 \item Die Kondition ist deshalb ($A$ ist symmetrisch)
476
  \begin{equation*}
477
   \kappa(A)
478
479
480
   =
   \frac{\lambda_\text{max}}{\lambda_\text{min}}
   =
481
482
483
484
485
486
487
   \frac{\lambda_{n-1,n-1}}{\lambda_{1,1}}
   =
   \frac{\sin^2 \Big(\frac{1}{2} (n-1) \pi h \Big)}{\sin^2 \Big(\frac{1}{2} \pi h \Big)}
   =
   \frac{\cos^2 \big(\frac{1}{2} \pi h \big)}{\sin^2 \big(\frac{1}{2} \pi h \big)}
   =
   \frac{4}{(\pi h)^2} \big(1+\mathcal{O} (h^2) \big).
488
489
490
491
492
493
494
495
496
497
498
499
500
  \end{equation*}

 \item Die Kondition hängt außerdem ab von:
  \begin{itemize}
   \item der Qualität der Elemente~\cite{shewchuk:2002}

   \item Bei einer Gleichung mit ortsabhängigen Koeffizienten wie
    \begin{equation*}
     - \div (D(x) \grad u) = f
    \end{equation*}
    hängt sie auch von den Eigenwerten von $D$ ab, und wie unterschiedlich
    diese auf $\Omega$ werden.
  \end{itemize}
501
502
\end{itemize}

503
504
\emph{Konsequenz:} Wir brauchen Lösungsverfahren, deren Geschwindigkeit
möglichst unabhängig von $\kappa(A)$ ist.
505

506
507
508
509
510

\begin{exercise}
 Wie hängt die Kondition bei einem dreidimensionalen Gitter von $h$ ab?
\end{exercise}

511
512

\subsubsection{Dünnbesetztheit}
513

514
Mehr zu dünnen Matrizen findet man z.B.\ bei \citet{pissanetzky:1984}.
515

516
\begin{itemize}
517
\item Ein Matrixeintrag $A_{ij} \colonequals a(\varphi_i, \varphi_j)$ ist genau dann
518
519
520
521
522
523
524
525
526
527
528
  nicht Null, wenn sich die Träger von $\varphi_i$ und $\varphi_j$ überlappen.

\item D.h.\ der überwiegende Teil der Matrixeinträge ist Null.

\item \glqq Eine Matrix heißt dünnbesetzt wenn sie so viele Nullen enthält dass es sich lohnt,
   dies auszunutzen.\grqq

\item Die Anzahl der Nicht-Null-Einträge pro Zeile ist durch eine (kleine) Konstante beschränkt,
  die nicht von der Gitterauflösung abhängt.

\item Die Matrix enthält also nur $O(n)$ Einträge
529
530
531
532
533
534
535
536
537
538
539
540
541
542
\end{itemize}

\bigskip

\emph{Konsequenz}: Wir brauchen Algorithmen und Datenstrukturen, die die Dünnbesetztheit ausnutzen.

\medskip Beispiel: Matrix--Vektor-Multiplikation $v = Aw$.

  \begin{algorithm}[H]
   \SetAlgoLined
   % \caption{...}
   \For{alle Zeilen $i$}
   {
     $v_i = 0$\\
543
     \For{alle Spalten $j$ in denen $A_{ij} \neq 0$}
544
545
546
547
548
     {
       $v_i = v_i + A_{ij} w_j$
     }
   }
  \end{algorithm}
549
Das braucht nur $O(\# \text{Nichtnulleinträge})$ Operationen (statt $O(n^2)$).
550

551
552
\bigskip

553
Direkte Verfahren wie Gauß-Elimination oder Cholesky-Zerlegung sind i.A. zu teuer.
554

555
556
557
558
559
560
561
\begin{itemize}
 \item \emph{Erinnerung:} Gauß-Elimination braucht $O(n^3)$ Rechenoperationen!

 \item Wendet man das Gauß-Verfahren auf solch eine Matrix an, so entstehen bei den Zwischenschritten in der Matrix eine beträchtliche Anzahl von zusätzlichen Einträgen ("`fill-in"').

 \item Das Gauß-Verfahren ist deshalb nicht nur zu langsam, es braucht auch zu viel Speicher.
\end{itemize}
562

563
\subsection{Direkte Verfahren für dünnbesetzte Gleichungen}
564

565
566
\begin{itemize}
 \item Es gibt direkte Verfahren zum Lösen von dünnbesetzten Gleichungssystemen.
567

568
 \item Sehr kompliziert / sehr ausgefeilt.
569

570
 \item Die Geschwindigkeit solcher Löser hängt nicht von der Kondition der Matrix ab.
571

572
573
574
575
576
577
578
579
580
581
582
583
 \item Manchmal gibt es Probleme mit Rundungsfehler, wenn Matrizen sehr schlecht
   konditioniert sind.

 \item Für kleine Probleme (In 2021: ca.\ 50\,000 Unbekannte oder weniger)
   sind sie konkurrenzlos schnell.

 \item Danach steigt die Rechenzeit schnell an.

 \item Limitierender Faktor ist aber häufig der Speicherverbrauch.

 \item Dennoch in der Praxis meist die Methode der Wahl.
\end{itemize}
584
585
586
587


\subsection{Lineare iterative Verfahren}

588
589
590
Iterative Verfahren konstruieren eine Folge $\{ x^k \}_{k \in \N}$,
die (hoffentlich) gegen die Lösung $x^*$ von $Ax=b$ konvergiert.

591
\bigskip
592

593
594
Für \emph{lineare} iterative Verfahren~\cite{dahmen_reusken:2008} schreibt man
$Ax=b$ als Fixpunktgleichung
595
\begin{equation*}
596
 x=x+C(b-Ax)
597
\end{equation*}
598
mit einer nichtsingulären Matrix $C \in \R^{n \times n}$.
599

600
\medskip
601

602
Die Größe $r(x) \colonequals b - Ax$
603
nennt man \emph{Residuum}.  Die Matrix $C$ heißt \emph{Vorkonditionierer}.
604

605
606
\bigskip

607
608
Dafür machen wir jetzt eine Fixpunktiteration:
\begin{equation*}
609
  x^{k+1}\colonequals x^k+C (b-Ax^k)=(I-CA)x^k+Cb,
610
611
612
613
  \qquad
  k = 0, 1, 2, \dots
\end{equation*}

614
Der Fehler im $k$-ten Iterationsschritt ist $e^k \colonequals x^k-x^*$. Es gilt
615
616
617
\begin{equation*}
 e^{k+1}=x^{k+1}-x^*
 =
618
\underbrace{(I-CA)}_{\text{Iterationsmatrix}}e^k.
619
\end{equation*}
620
Das erklärt den Namen: Die Abbildung $e^k \mapsto e^{k+1}$ ist linear.
621
622
623
624
625

\begin{definition}
 Die Matrix $I - CA$ heißt \emph{Iterationsmatrix} der Methode.
\end{definition}

626
\subsubsection{Konvergenz:}
627
628

\begin{theorem}
629
 Sei $\norm{\cdot}$ eine submultiplikative Matrixnorm.
630
631
 Das Verfahren konvergiert für jeden Startwert $x^0 \in \R^n$ gegen die Lösung
 von $Ax = b$, wenn $\norm{I - CA} < 1$.
632
633
\end{theorem}

634
635
636
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
637
638
639
640
641
642
643
644
645
\begin{equation*}
 \abs{\lambda}\norm{v}
 =
 \norm{\lambda v}
 =
 \norm{B v}
 \le
 \norm{B}\norm{v}.
\end{equation*}
646
647

Damit erhält man:
648

649
650
651
\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}
652

653
654
655
656
657
658
659
660
661
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:
662
\begin{itemize}
663
 \item Das Jacobi-Verfahren
664
  \begin{equation*}
665
666
   C = D^{-1}
  \end{equation*}
667

668
669
670
671
 \item Das Gauß--Seidel-Verfahren
  \begin{equation*}
   C = (D+L)^{-1}
  \end{equation*}
672

673
 \item Das symmetrische Gauß--Seidel-Verfahren
674
  \begin{equation*}
675
   C = (D+L)^{-1}+(D+R)^{-1}-(D+R)^{-1}A(D+L)^{-1}
676
677
678
679
680
681
682
  \end{equation*}
\end{itemize}

\emph{Beachte:} Die Matrix $C$ wird nie explizit ausgerechnet!


\subsection{Das Jacobi-Verfahren}
683
684
685
686

Wir betrachten dieses Verfahren etwas genauer.  Man kann es auch
\glqq direkt\grqq motivieren.

687
Wir nehmen im Folgenden an, dass $A_{ii} \neq 0$ für alle $i=1,\ldots,n$.
688
689
690
691
692

\medskip

Betrachte das lineare Gleichungssystem:
\begin{align*}
693
694
695
    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
696
\end{align*}
697
\emph{Idee:} Löse $i$-te Zeile nach $x_i$ auf für $i=1,\ldots,n$:
698
\begin{align*}
699
700
701
    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
702
703
704
\end{align*}
Mache daraus ein iteratives Verfahren:
\begin{align*}
705
706
707
    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
708
709
710
\end{align*}
Für $i=1,\ldots,n$
\begin{equation*}
711
712
 x_i^{k+1}
 =
713
 \frac{1}{A_{ii}} \bigg(b_i-\sum_{\substack{j=1\\j\neq i}}^n A_{ij}x_j^k \bigg)
714
 =
715
 x^k_i + \frac{1}{A_{ii}} \bigg(b_i-\sum_{j=1}^n A_{ij}x_j^k \bigg).
716
\end{equation*}
717
718
719
Beachte:
\begin{itemize}
 \item Die Rechnungen für die verschiedenen $i$ sind voneinander unabhängig $\implies$ leicht zu parallelisieren.
720

721
722
 \item In der Praxis gilt die Summe natürlich nur über die Einträge der $i$-ten Zeile von $A$, die $\neq 0$ sind.
\end{itemize}
723

724
725
726
727
\begin{exercise}
 Rechnen Sie nach dass dies wirklich ein lineares Verfahren mit Vorkonditionierer
 $C = D^{-1}$ ist.
\end{exercise}
728

729
730
731
732
733
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}
 =
734
 x^k_i + \frac{\eta}{A_{ii}} \bigg(b_i-\sum_{j=1}^n A_{ij}x_j^k \bigg)
735
736
737
 \qquad \text{bzw.} \qquad
 x^{k+1} = x^k + \eta D^{-1} ( b - Ax^k).
\end{equation*}
738
739
740
741
742

\subsubsection{Konvergenz}

Wir wenden das bekannte Konvergenzkriterium an:
\begin{theorem}
743
Das Jacobi Verfahren konvergiert genau dann, wenn $\rho (I- \eta D^{-1}A ) <1$.
744
745
\end{theorem}

746
747
748
749
750
751
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}.
752

753
754
Wir greifen die Konvergenztheorie des Jacobi-Verfahrens in Kapitel~\ref{sec:jacobi_smoother}
wieder auf.
755

756
\subsubsection{Konvergenzgeschwindigkeit}
757
758
759
760
761

Die folgende Rechnung stammt aus~\citet[Beispiel~13.10]{dahmen_reusken:2008}.

\bigskip

762
763
Sei $A$ die Matrix des Poisson-Problems auf einem uniformen Gitter
mit Dreiecks- oder Viereckselementen.  Eine Zeile von $A$ ist dann
764
\begin{equation*}
765
    \frac{1}{h^2} (4u_{i,j}-u_{i-1,j}-u_{i+1,j}-u_{i,j-1}-u_{i,j+1} )=f_i.
766
767
\end{equation*}

768
769
770
771
\begin{exercise}
 Rechnen Sie nach dass man tatsächlich diese Matrix erhält.
\end{exercise}

772

773
774
775

Jacobi-Verfahren: $D=4h^{-2}I$.

776
Wir versuchen jetzt, den Spektralradius von $I-CA = I - \eta D^{-1}A$ abzuschätzen:
777
\begin{align*}
778
779
780
781
    \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}}.
782
783
\end{align*}
Der Bruch ist gerade der \emph{Rayleigh-Quotient}
784
785
786
787
788
789
\begin{equation*}
  R(A,x) \colonequals \frac{x^TAx}{\Vert x \Vert^2}.
\end{equation*}
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.

790
791
\medskip

792
Daraus folgt dass
793
794
795
\begin{align*}
 \rho (I - \eta D^{-1}A )
 & =
796
 \sup \lbrace \abs[\Big]{1- \frac{1}{4} \eta h^2 \lambda} \; : \; \text{$\lambda$ Eigenwert von $A$} \rbrace \\
797
798
 %
 & =
799
 1-2 \eta \sin^2 ( \frac{1}{2} \pi h ).
800
\end{align*}
801
802
Taylor-Entwicklung:
\begin{equation*}
803
    \sin^2 s = s^2 - \frac{1}{3}s^4 + \frac{2}{45}s^6 + \cdots
804
\end{equation*}
805
ergibt
806
\begin{equation*}
807
 \frac{\norm{e^{k+1}}}{\norm{e^k}} \approx \rho (I-D^{-1}A) \approx 1-\frac{\eta}{2} \pi^2h^2
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
\end{equation*}
Dieser Ausdruck geht "`quadratisch"' (also ziemlich schnell) gegen $1$ wenn $h \to 0$.

\bigskip

Wir wollen ausrechnen, wie viele Iterationen man ungefähr braucht, um den Anfangsfehler
um einen Faktor $R$ zu reduzieren.  D.h., welches $k$ soll man wählen, um ungefähr
\begin{equation*}
 \frac{\norm{e^k}}{\norm{e^0}} \le \frac{1}{R}
\end{equation*}
zu erhalten?

Da $\frac{\norm{e^k}}{\norm{e^0}} \approx \rho^k$ erhält man
\begin{equation*}
 k
 =
 \log_\rho \frac{1}{R}
 =
826
 \frac{-\ln R}{\ln \rho (I-D^{-1}A )} \approx \frac{-\ln R}{\ln (1-\frac{1}{2}\pi^2 h^2 )} \approx \frac{2}{\pi^2 h^2} \ln R,
827
828
829
\end{equation*}
da $\ln x \approx (x-1)-\frac{1}{2} (x-1)^2 + \ldots$ ist.

830
Für ein doppelt so feines Gitter braucht man also viermal so viele Iterationen
831
832
(und diese sind natürlich auch noch teurer.).

833
834
\medskip

835
Das Jacobi-Verfahren scheint also kein sonderlich gutes Verfahren zu sein.
836
837
838
839
840
841
842

\subsection{Das Gauß-Seidel-Verfahren}

Betrachte noch einmal die Jacobi-Rechenvorschrift:
\begin{equation*}
 x_i^{k+1}
 =
843
844
 \frac{1}{A_{ii}} (b_i-A_{i1}x_1^{k} - A_{i2}x_2^{k}-
   \ldots -A_{i,i-1}x_{i-1}^k - A_{i,i+1}x_{i+1}^k - \ldots - A_{in}x_n^k ).
845
846
847
\end{equation*}
Eigentlich haben wir für $x_1,\ldots,x_{i-1}$ schon bessere Werte als $x_1^k,\ldots,x_{i-1}^k$, nämlich $x_1^{k+1},\ldots,x_{i-1}^{k+1}$.

Sander, Oliver's avatar
Sander, Oliver committed
848
849
\medskip

850
851
Gauß-Seidel-Verfahren:
\begin{align*}
Sander, Oliver's avatar
Sander, Oliver committed
852
853
 x_i^{k+1}
 & =
854
 \frac{1}{A_{ii}} (b_i-A_{i1}x_1^{k+1} - \ldots -A_{i,i-1}x_{i-1}^{k+1}-A_{i,i+1}x_{i+1}^k-\ldots-A_{in}x_n^k ) \\
Sander, Oliver's avatar
Sander, Oliver committed
855
 & =
856
 \frac{1}{A_{ii}} (b_i-\sum_{j=1}^{i-1} A_{ij}x_j^{k+1} -\sum_{j=i+1}^n A_{ij} x_j^k ).
857
858
859
\end{align*}
$x_i^{k+1}$ hängt von $x_{i-1}^{k+1}$ ab $\implies$ keine Parallelisierung möglich.

860
861
862
863
864
\begin{exercise}
 Seien $D$, $L$ Diagonalteil bzw.\ linker Dreiecksteil von $A$.
 Rechnen Sie nach dass das Gauß-Seidel-Verfahren ein lineares Verfahren
 mit $C = (D+L)^{-1}$ ist.
\end{exercise}
865
866
867
868
869
870
871
872
873

\subsubsection{Konvergenz}

Wir erwarten bessere Konvergenzeigenschaften als für das Jacobi-Verfahren.

\medskip

Und in der Tat:
\begin{theorem}
874
875
876
877
Das Gauß-Seidel-Verfahren konvergiert, wenn
\begin{itemize}
 \item $A$ symmetrisch und positiv definit ist und/oder
 \item $A$ irreduzibel und diagonaldominant ist.
878
879
880
881
882
883
\end{itemize}
\end{theorem}
\subsubsection{Konvergenzgeschwindigkeit}
\begin{example}
Sei $A$ die Matrix des Poisson-Problems für ein quadratisches Gebiet.
\begin{itemize}
Sander, Oliver's avatar
Sander, Oliver committed
884
 \item Es gilt $\rho(I-(D+L)^{-1}A) = ( \rho (I-D^{-1}A ) )^2$.
885
886
887
888
  Das heißt der Spektralradius der Jacobi-Methode ist das Quadrat
  des Spektralradius der Gauß-Seidel Methode.  Laut \citet{dahmen_reusken:2008}
  steht das bei \cite{hackbusch:1994}.

889
890
  \todo[inline]{I overlooked Remark 5.5.2. in [3] (On the top of page 131). There is the explicit statement.}

891
892
893
894
895
  \begin{exercise}
   Finden Sie die genaue Stelle, und sagen Sie mir bescheid.
  \end{exercise}


896
897
898
 %
 \item Deshalb
  \begin{equation*}
Sander, Oliver's avatar
Sander, Oliver committed
899
   \rho(I-(D+L)^{-1}A) = \cos^2(\pi h)
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
  \end{equation*}

 \item Taylor-Entwicklung für $\cos^2$
  \begin{equation*}
   \cos^2 \pi h
   \approx
   \Big(1-\frac{1}{2} \pi^2h^2+\ldots \Big)^2
   =
   1-2 \frac{1}{2} \pi^2h^2+\frac{1}{4} \pi^4h^4+\ldots \approx 1-\pi^2h^2
  \end{equation*}
 %
 \item Fehlerreduktion
 \begin{equation*}
  \frac{\norm{e^{k+1}}}{\norm{e^k}} \approx 1-\pi^2h^2
 \end{equation*}
 %
 \item Um den Startfehler um den Faktor $R$ zu reduzieren, braucht man etwa
 \begin{equation*}
  \frac{- \ln R}{\ln \rho (I-(D-L)^{-1}A)}
   \approx
920
  \frac{- \ln R}{\ln (1-\pi^2h^2 )} \approx \frac{\ln R}{\pi^2h^2}
921
922
923
924
925
926
 \end{equation*}
 Iterationen.
\end{itemize}
\end{example}
\emph{Faustregel:} Das Gauß-Seidel-Verfahren braucht nur etwa halb so viele Iterationen wie das Jacobi-Verfahren.

927
928
929
930
931
932
933
934
935
\subsection{Nichtlineare Verfahren}

Zum Vergleich sind hier noch Konvergenzresultate für zwei bekannte
nichtlineare Verfahren.  (Beide ausschließlich für symmetrische Matrizen).

Die Konvergenzresultate sind als Funktion der Konditionszahl $\kappa$
formuliert.

\subsubsection{Das Gradientenverfahren}
936
937
938
939
940
941
942
943
944

\begin{theorem}[{\cite[Kapitel~6]{shewchuk:1994}}]
Sei $\kappa$ die Kondition der Matrix $A$. Dann gilt nach $k$ Schritten des Gradientenverfahrens
\begin{equation*}
 \norm{e^k}_A \leq \Big(\frac{\kappa-1}{\kappa+1}\Big)^k \norm{e^0}_A,
\end{equation*}
wobei $\norm{\cdot}_A$ die Energienorm ist.
\end{theorem}

945
\subsubsection{Das Verfahren der konjugierten Gradienten (CG-Verfahren)}
946
947

\begin{itemize}
948
949
950
951
952
953
 \item Zuerst vorgeschlagen von \citet{hestenes_stiefel:1952} im Jahr 1952.
   Eine gute Beschreibung findet sich bei~\citet{shewchuk:1994}.
 \item Toll: ein direkter Löser für dünnbesetzte lineare GS mit Komplexität $O(n \cdot$ Anzahl der Einträge von $A)$.
 \item Funktioniert in der Praxis aber schlecht: Wegen Rundungsfehlern erhält man
    häufig nach $n$ Iterationen \emph{nicht} die Lösung.

954
955
956
957
 \item Geriet zwischenzeitlich in Vergessenheit
 \item Erlebte Revival als iterative Methode
\end{itemize}

958
959
\begin{theorem}[{\cite{shewchuk:1994}}]
Sei $\kappa$ die Kondition der Matrix $A$. Dann gilt nach $k$ Schritten des CG-Verfahrens
960
\begin{align*}
961
 \norm{e^k}_A \leq 2 ( \frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1} )^k \cdot \norm{e^0}_A.
962
\end{align*}
963
964
\end{theorem}

965
966
967
\begin{itemize}
 \item Diese Abschätzung für CG ist aber sehr schwach. In vielen Fällen konvergiert das Verfahren deutlich besser!

968
 \item Das CG-Verfahren wird tatsächlich sehr häufig für FE-Probleme verwendet.
969

970
971
 \item Aber man benötigt wieder einen Vorkonditionierer, sonst ist auch
   das CG-Verfahren zu langsam.
972

973
 \item Dabei wird wieder die Gleichung $0 = b - Ax$ durch $0 = C(b-Ax)$ ersetzt.
974

975
976
 \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.
977

978
979
 \item Die Konstruktion passender $C$ bildet ein großes Gebiet innerhalb
   der Numerik.
980

981
 \item Eine Möglichkeit: Mehrgitter!
982
983
\end{itemize}

984
\section{Mehrgitterverfahren: Der anschauliche Zugang}
985

986
987
988
Das Ziel ist es, einen iterativen Löser zu konstruieren, dessen Konvergenzgeschwindigkeit
nicht von der Gitterweite $h$ abhängt.

989
\subsection{Glättung}
990

991
Gauß--Seidel- und Jacobi-Verfahren sind für typische Finite-Elemente-Matrizen
992
extrem langsam.
993

994
995
996
997
\medskip

Sie haben allerdings eine wichtige Eigenschaft, die wir uns
für Mehrgitterverfahren zunutze machen:
998

999
1000
1001
\smallskip

$\longrightarrow$ Die zwei Verfahren \emph{glätten}.
1002

1003
1004
1005
\bigskip

\emph{Beispiel:} Wir betrachten wieder das Poisson-Problem
1006
1007
1008
1009
1010
\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$.
1011
1012
1013

\medskip

1014
1015
1016
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:
1017
1018
1019
1020
1021
1022
1023
\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.

1024
\medskip
1025
1026
1027

Als Startiterierte $x^0$ wählen wir $x^0 = 0 \in \R^n$.

1028
\medskip
1029
1030
1031
1032

Hier sind einige der ersten Iterierten:

\begin{center}
1033
 \begin{tabular}{ccc}
1034
1035

  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_zero_iterate0} &
1036
1037
  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_zero_iterate1} &
  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_zero_iterate2} \\
1038

1039
  $x^0$ & $x^1$ & $x^2$ \\
1040

1041
  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_zero_iterate3} &
1042
1043
1044
  \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} \\

1045
  $x^3$ & $x^{10}$ & $x^{50}$
1046
1047
1048
1049
1050
1051
1052
1053
 \end{tabular}
\end{center}

Man erkennt dass das Verfahren gehen die Lösung konvergiert, aber es braucht
tatsächlich recht lange.

\bigskip

1054
Jetzt zeigen wir was wir mit \glqq Gauß--Seidel glättet\grqq{} meinen.
1055
1056
1057

\medskip

1058
1059
1060
1061
Wir lösen nochmal die Gleichung, fangen jetzt aber von einem
verrauschten Startwert an.

\begin{center}
1062
 \begin{tabular}{ccc}
1063
  \includegraphics[width=0.3\textwidth,trim=0 50 0 0]{gs_from_noise_iterate0} &
1064
1065
  \includegraphics[width=0.3\textwidth,trim=0 50 0 0]{gs_from_noise_iterate1} &
  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_noise_iterate2} \\
1066

1067
  $x^0$ & $x^1$ & $x^2$ \\
1068

1069
  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_noise_iterate3} &
1070
1071
1072
  \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} \\

1073
  $x^3$ & $x^{10}$ & $x^{50}$
1074
1075
1076
 \end{tabular}
\end{center}

1077
\begin{itemize}
1078
 \item Das Verfahren konvergiert immer noch.
1079

1080
 \item Es scheint weder schneller noch langsamer geworden zu sein.
1081

1082
1083
1084
 \item Das Rauschen ist aber schon nach sehr wenigen Iterationen weitestgehend
   verschwunden.
\end{itemize}
1085

1086
1087
1088
1089
1090
1091
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}}
Sander, Oliver's avatar
Sander, Oliver committed
1092
   \bigg( b_i - \sum_{j=1}^{i-1} A_{ij} x^{k+1}_j - \sum_{j=i+1}^{n} A_{ij} x^k_j \bigg).
1093
1094
1095
1096
\end{equation*}
Dabei tauchen in den Summen nur die Nachbarknoten auf.
\begin{itemize}
 \item Zumindest in diesem Modellproblem sind die entsprechenden Koeffizienten
1097
   (fast) alle gleich.
1098
1099
1100
1101
1102
1103

 \item $x^k_i$ wird also durch den Mittelwert seiner Nachbarn ersetzt.
\end{itemize}

\bigskip

1104
So richtig überzeugend ist diese Argumentation trotzdem noch nicht:
1105

1106
1107
1108
\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?
1109

1110
 \item In so einem Fall könnte doch das Gauß--Seidel-Verfahren nicht glätten?
1111
1112
\end{itemize}

1113
1114
In Wirklichkeit ist es etwas anders:

1115
\begin{itemize}
1116
1117
1118
 \item Das Gauß--Seidel-Verfahren glättet \emph{nicht} die Iterierten.

 \item Das Gauß--Seidel-Verfahren glättet die \emph{Fehler}!
1119
1120
\end{itemize}

1121
\missingfigure{Sechs Bilder mit den Fehlern für die obigen Iterierten!}
1122

1123
1124
1125
\subsection{Fourier-Analyse}
\label{sec:jacobi_smoother}

1126
1127
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.
1128

1129
1130
\medskip

1131
Wir untersuchen hier das Jacobi-Verfahren.  Für das Gauß--Seidel-Verfahren geht es auch,
1132
ist aber etwas komplizierter.  Interessierte finden die Rechnung zum Beispiel bei
1133
1134
1135
\citet[Kapitel~4]{trottenberg_oosterlee_schueller:2001}

\bigskip
1136
1137
1138
1139
1140

Zur Erinnerung: Das gedämpfte Jacobi-Verfahren ist
\begin{equation*}
 x^{k+1}
 =
1141
1142
1143
 x^k + \eta D^{-1} (b - Ax^k)
 =
 (I - \eta D^{-1}A)x^k + \eta D^{-1}b.
1144
1145
1146
1147
\end{equation*}

Wir betrachten wieder das Modellproblem.  Dort ist
\begin{equation*}
1148
 D_i \colonequals A_{ii} = \frac{4}{h^2}
1149
 \qquad
1150
 \text{für alle $i=1,\dots,N$},
1151
1152
1153
\end{equation*}
und somit
\begin{equation*}
1154
 x^{k+1} = \Big(I- \eta \frac{h^2}{4} A \Big) x^k + \eta \frac{h^2}{4}b.
1155
1156
1157
1158
1159
\end{equation*}

Um das Glättungsverhalten des Jacobi-Verfahrens zu verstehen berechnen wir
die Eigenfunktionen der Iterationsmatrix $I - \omega \frac{h^2}{4}A$.

1160
1161
\medskip

1162
1163
1164
Die Eigenfunktionen von $I - \omega \frac{h^2}{4}A$ sind offensichtlich
die selben wie die von $A$, und zwar
\begin{equation*}
1165
 \varphi_h^{\nu,\mu}
1166
 =
1167
1168
1169
 \sin \nu\pi x \cdot \sin \mu \pi y
 \qquad
 1 \le \nu,\mu < n \colonequals \frac{1}{h}.
1170
1171
1172
\end{equation*}
Die Eigenwerte von $A$ sind
\begin{equation*}
1173
 \lambda^{\nu\mu}
1174
 =