skript-mehrgitter-sander.tex 323 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
\usetikzlibrary{calc}
23
\usepackage{pgfplots}
Sander, Oliver's avatar
Sander, Oliver committed
24
\pgfplotsset{compat=1.16}
25
26
27
28

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

39
\usepackage[maxbibnames=99,  % Never write 'et al.' in a bibliography
40
41
42
43
44
45
46
            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}
47
\usepackage[disable]{todonotes}
48
49
50
51
52
53
54
55
\usepackage{esdiff}
\usepackage{hyperref}

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

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

65
66
67
\newcommand{\Radd}{R^\textup{add}}
\newcommand{\Rmult}{R^\textup{mult}}

68
% Bold letters
69
\newcommand*{\ba}{\bm{a}}
70
71
72
\newcommand*{\bb}{\bm{b}}
\newcommand*{\bd}{\bm{d}}
\newcommand*{\be}{\bm{e}}
73
\newcommand*{\bof}{\bm{f}}
74
\newcommand*{\bh}{\bm{h}}
75
76
\newcommand*{\bj}{\bm{j}}
\newcommand*{\bn}{\mathbf{n}}
77
78
79
\newcommand*{\bp}{\bm{p}}
\newcommand*{\bq}{\bm{q}}
\newcommand*{\bu}{\bm{u}}
80
\newcommand*{\bv}{\mathbf{v}}
81
82
\newcommand*{\bx}{\bm{x}}
\newcommand*{\by}{\bm{y}}
83
\newcommand*{\boz}{\bm{z}}
84

85

86
\newcommand*\bA{\bm{A}}
87
88
89
90
\newcommand*\bB{\bm{B}}
\newcommand*\bD{\bm{D}}
\newcommand*\bE{\bm{E}}
\newcommand*\bF{\bm{F}}
91
92
93
94
95
96
\newcommand*\bH{\bm{H}}
\newcommand*\bI{\bm{I}}
\newcommand*\bL{\bm{L}}
\newcommand*\bP{\bm{P}}
\newcommand*\bQ{\bm{Q}}
\newcommand*\bV{\bm{V}}
97
\newcommand*\bZ{\bm{Z}}
98

99
\newcommand*\bpsi{\bm{\psi}}
100
101
\newcommand{\btau}{\bm{\tau}}

102
103
\newcommand*{\bPi}{\bm{\Pi}}

104
105
106
107
108
109
110
111
112
113
114
115
\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}}
116
117
\DeclarePairedDelimiter{\abs}{\lvert}{\rvert}
\DeclarePairedDelimiter{\norm}{\lVert}{\rVert}
118
\DeclarePairedDelimiter{\triplenorm}{\vvvert}{\vvvert}
119
120
121

\theoremstyle{plain} %Text ist Kursiv
\newtheorem{theorem}{Satz}[chapter]
122
\newtheorem{hilfssatz}[theorem]{Hilfssatz}
123
124
125
126
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Korollar}
\newtheorem{definition}[theorem]{Definition}
127
\newtheorem*{definition*}{Definition}
128
129
130
131
132
\newtheorem{problem}[theorem]{Problem}

\theoremstyle{definition} %Text ist \"upright"
\newtheorem{remark}[theorem]{Bemerkung}
\newtheorem{example}[theorem]{Beispiel}
133
\newtheorem{exercise}[theorem]{Übung}
134

135
136
137
138
% Silbentrennung für Worte die TeX nicht alleine kann
% TODO: Reicht es eventuell zu sagen dass dies ein englisches Word ist?
\hyphenation{strength-ened}

139
140
141
142
143
144
145
146
147
148
149
150
151
\swapnumbers

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

\begin{document}

\begin{titlepage}

\begin{center}

\vspace*{0.3\textheight}

152
{\Huge Mehrgitterverfahren}
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
\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
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
\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$.
211
Wir betrachten als Beispiel die \emph{Poisson-Gleichung}
Sander, Oliver's avatar
Sander, Oliver committed
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
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
\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}
262
 a(u,v) = \ell(v)
Sander, Oliver's avatar
Sander, Oliver committed
263
264
265
\end{equation}
mit
\begin{equation*}
266
267
268
 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
269
270
271
+
\int_\Omega f v \, dx.
\end{equation*}
272
273
Man nennt \eqref{eq:schwache_formulierung} die \emph{schwache Formulierung}
(oder \emph{Variationsformulierung})
Sander, Oliver's avatar
Sander, Oliver committed
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
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
295
a(u,v) = \ell(v)
Sander, Oliver's avatar
Sander, Oliver committed
296
297
\qquad v \in H
\end{equation*}
298
for jedes $\ell \in H'$ eine eindeutig bestimmte Lösung.
Sander, Oliver's avatar
Sander, Oliver committed
299
300
301
302
303
304
305
306
307
\end{theorem}




\section{Diskretisierung (Galerkin-Verfahren)}

Sei $V_h$ ein endlichdimensionaler Teilraum von $H^1_{\Gammatight{D}}$ (mit Dimension $n$).
Der Index $h$
308
bezeichne einen Ordnungsparameter.  Die \emph{diskrete Variationsformulierung}
Sander, Oliver's avatar
Sander, Oliver committed
309
310
311
lautet
\begin{equation}
\label{eq:diskrete_formulierung}
312
 u_h \in V_h \quad : \quad a(u_h,v_h) = \ell(v_h)
Sander, Oliver's avatar
Sander, Oliver committed
313
314
315
316
317
318
319
320
321
\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}
322
 u_h \in V_h \quad : \quad a(u_h, \varphi_i) = \ell(\varphi_i)
Sander, Oliver's avatar
Sander, Oliver committed
323
324
325
\qquad \text{für alle $0 \le i < n$}.
\end{equation*}
Dies wiederum entspricht dem linearen Gleichungssystem
326
327
\begin{equation}
\label{eq:fe_lineares_gleichungssystem}
328
 Ax = b,
329
\end{equation}
Sander, Oliver's avatar
Sander, Oliver committed
330
331
332
wobei
\begin{alignat*}{2}
A &\in \R^{n \times n},
333
\qquad &A_{ij} & \colonequals a(\varphi_i, \varphi_j)
Sander, Oliver's avatar
Sander, Oliver committed
334
335
336
337
338
=
\int_\Omega \nabla \varphi_i \nabla \varphi_j \, dx \\
%
b &\in \R^n,
\qquad &
339
b_i & \colonequals \ell(\varphi_i)
Sander, Oliver's avatar
Sander, Oliver committed
340
341
342
343
344
=
\int_{\Gamma_N} g \varphi_i \, ds
+
\int_\Omega f \varphi_i \, dx
\end{alignat*}
345
und $x \in \R^n$ die Koeffizienten von $u_h$ bzgl.\ der Basis $\{\varphi_i\}$ sind.
Sander, Oliver's avatar
Sander, Oliver committed
346

347
348
349
\bigskip

Die Matrix $A$ nennt man \emph{Steifigkeitsmatrix}.  Dieser Ausdruck kommt aus den
Sander, Oliver's avatar
Sander, Oliver committed
350
351
352
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$
353
nennt man deshalb manchmal auch \emph{Lastvektor}.
Sander, Oliver's avatar
Sander, Oliver committed
354
355
356
357
358
359
360
361


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

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

362
363
364
365
\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
366
367
368
369
370
371
372
373
374
375
376
377
378
379
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}$
380
von abgeschlossenen Dreiecken $T$ heißt \emph{Triangulierung} von $\Omega$, falls gilt
Sander, Oliver's avatar
Sander, Oliver committed
381
382
383
\begin{enumerate}
 \item Die Menge der Dreiecke überdeckt den Abschluss des Gebiets
\begin{equation*}
384
 \overline{\Omega} = \bigcup_{T \in \mathcal{T}} T.
Sander, Oliver's avatar
Sander, Oliver committed
385
386
387
388
389
390
391
392
\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*}
393
 h \colonequals \max_{T \in \mathcal{T}} \operatorname{diam} T,
Sander, Oliver's avatar
Sander, Oliver committed
394
395
\end{equation*}
der als Index in $V_h$ auftaucht, beschreibt traditionell die Feinheit des Gitters.
396
Oft nennt man $h$ auch \emph{Gitterweite}.
Sander, Oliver's avatar
Sander, Oliver committed
397
398
399
400

\begin{definition}
Der Raum
\begin{equation*}
401
 V_h^{p}
Sander, Oliver's avatar
Sander, Oliver committed
402
403
404
\colonequals
\{ v \in C(\overline{\Omega})
\; | \;
405
\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
406
\end{equation*}
407
heißt Raum der Lagrange-Elemente $p$-Ordnung
Sander, Oliver's avatar
Sander, Oliver committed
408
409
410
411
412
413
414
415
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}
416
\caption{Links: eine Finite-Elemente Funktion.  Rechts: ein Element der Knotenbasis erster Ordnung}
Sander, Oliver's avatar
Sander, Oliver committed
417
418
419
\label{fig:beispiel-fe-funktion}
\end{figure}

420
421
422
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
423
424
425
\begin{equation*}
 \varphi_i (v_j) = \delta_{ij}
\qquad
426
\text{für alle $v_j \in \mathcal{L}$}.
Sander, Oliver's avatar
Sander, Oliver committed
427
\end{equation*}
428
Es ist also insbesondere $\dim V_h^{p} = \abs{\mathcal L \setminus \mathcal{L}_{\Gamma_D}}$
Sander, Oliver's avatar
Sander, Oliver committed
429
430
(Abbildung~\ref{fig:beispiel-fe-funktion}).

431
432
433
434
\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
435
436

\begin{theorem}
437
438
439
 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
440
\begin{equation*}
441
 \norm{u - u_h}_1 \le C h^p \abs{u}_{p+1}
Sander, Oliver's avatar
Sander, Oliver committed
442
443
444
445
446
447
448
449
\end{equation*}
mit einer Konstanten $C$.
\end{theorem}

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

450
\chapter{Mehrgitterverfahren}
Sander, Oliver's avatar
Sander, Oliver committed
451

452
\section{Iterative Lösungsverfahren für dünnbesetzte lineare Gleichungssysteme}
453
\label{sec:iterative_loeser}
454

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

458
\subsection{Eigenschaften von Finite-Elemente-Matrizen}
459
460
461

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

462
\subsubsection{Größe}
463
Finite-Elemente-Matrizen können sehr groß werden:
464
\begin{itemize}
465
466
467
468
469
470
471
 \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.

472
 \item Eine Zahl in doppelter Genauigkeit braucht $8$~Byte.  In ein Gigabyte Hauptspeicher
473
474
475
476
477
478
479
480
481
   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}

482
483
Die Matrizen zu elliptischen Problemen sind positiv definit, aber sie
sind schlecht konditioniert.
484

485
\medskip
486

487
488
489
490
491
492
493
494
495
496
497
498
499
\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)
500
  \begin{equation*}
501
   \kappa(A)
502
503
504
   =
   \frac{\lambda_\text{max}}{\lambda_\text{min}}
   =
505
506
507
508
509
510
511
   \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).
512
513
514
515
516
517
518
519
520
521
522
523
524
  \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}
525
526
\end{itemize}

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

530
531
532
533
534

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

535
536

\subsubsection{Dünnbesetztheit}
537

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

540
\begin{itemize}
541
\item Ein Matrixeintrag $A_{ij} \colonequals a(\varphi_i, \varphi_j)$ ist genau dann
542
543
544
545
546
547
548
549
550
551
552
  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
553
554
555
556
557
558
559
560
561
562
563
564
565
566
\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$\\
567
     \For{alle Spalten $j$ in denen $A_{ij} \neq 0$}
568
569
570
571
572
     {
       $v_i = v_i + A_{ij} w_j$
     }
   }
  \end{algorithm}
573
Das braucht nur $O(\# \text{Nichtnulleinträge})$ Operationen (statt $O(n^2)$).
574

575
576
\bigskip

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

579
580
581
582
583
584
585
\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}
586

587
\subsection{Direkte Verfahren für dünnbesetzte Gleichungen}
588

589
590
\begin{itemize}
 \item Es gibt direkte Verfahren zum Lösen von dünnbesetzten Gleichungssystemen.
591

592
 \item Sehr kompliziert / sehr ausgefeilt.
593

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

596
597
598
599
600
601
602
603
604
605
606
607
 \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}
608
609
610
611


\subsection{Lineare iterative Verfahren}

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

615
\bigskip
616

617
618
Für \emph{lineare} iterative Verfahren~\cite{dahmen_reusken:2008} schreibt man
$Ax=b$ als Fixpunktgleichung
619
\begin{equation*}
620
 x=x+C(b-Ax)
621
\end{equation*}
622
mit einer nichtsingulären Matrix $C \in \R^{n \times n}$.
623

624
\medskip
625

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

629
630
\bigskip

631
632
Dafür machen wir jetzt eine Fixpunktiteration:
\begin{equation*}
633
  x^{k+1}\colonequals x^k+C (b-Ax^k)=(I-CA)x^k+Cb,
634
635
636
637
  \qquad
  k = 0, 1, 2, \dots
\end{equation*}

638
Der Fehler im $k$-ten Iterationsschritt ist $e^k \colonequals x^k-x^*$. Es gilt
639
640
641
\begin{equation*}
 e^{k+1}=x^{k+1}-x^*
 =
642
\underbrace{(I-CA)}_{\text{Iterationsmatrix}}e^k.
643
\end{equation*}
644
Das erklärt den Namen: Die Abbildung $e^k \mapsto e^{k+1}$ ist linear.
645
646
647
648
649

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

650
\subsubsection{Konvergenz:}
651
652

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

658
659
660
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
661
662
663
664
665
666
667
668
669
\begin{equation*}
 \abs{\lambda}\norm{v}
 =
 \norm{\lambda v}
 =
 \norm{B v}
 \le
 \norm{B}\norm{v}.
\end{equation*}
670
671

Damit erhält man:
672

673
674
675
\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}
676

677
678
679
680
681
682
683
684
685
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:
686
\begin{itemize}
687
 \item Das Jacobi-Verfahren
688
  \begin{equation*}
689
690
   C = D^{-1}
  \end{equation*}
691

692
693
694
695
 \item Das Gauß--Seidel-Verfahren
  \begin{equation*}
   C = (D+L)^{-1}
  \end{equation*}
696

697
 \item Das symmetrische Gauß--Seidel-Verfahren
698
  \begin{equation*}
699
   C = (D+L)^{-1}+(D+R)^{-1}-(D+R)^{-1}A(D+L)^{-1}
700
701
702
703
704
705
706
  \end{equation*}
\end{itemize}

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


\subsection{Das Jacobi-Verfahren}
707
708
709
710

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

711
Wir nehmen im Folgenden an, dass $A_{ii} \neq 0$ für alle $i=1,\ldots,n$.
712
713
714
715
716

\medskip

Betrachte das lineare Gleichungssystem:
\begin{align*}
717
718
719
    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
720
\end{align*}
721
\emph{Idee:} Löse $i$-te Zeile nach $x_i$ auf für $i=1,\ldots,n$:
722
\begin{align*}
723
724
725
    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
726
727
728
\end{align*}
Mache daraus ein iteratives Verfahren:
\begin{align*}
729
730
731
    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
732
733
734
\end{align*}
Für $i=1,\ldots,n$
\begin{equation*}
735
736
 x_i^{k+1}
 =
737
 \frac{1}{A_{ii}} \bigg(b_i-\sum_{\substack{j=1\\j\neq i}}^n A_{ij}x_j^k \bigg)
738
 =
739
 x^k_i + \frac{1}{A_{ii}} \bigg(b_i-\sum_{j=1}^n A_{ij}x_j^k \bigg).
740
\end{equation*}
741
742
743
Beachte:
\begin{itemize}
 \item Die Rechnungen für die verschiedenen $i$ sind voneinander unabhängig $\implies$ leicht zu parallelisieren.
744

745
746
 \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}
747

748
749
750
751
\begin{exercise}
 Rechnen Sie nach dass dies wirklich ein lineares Verfahren mit Vorkonditionierer
 $C = D^{-1}$ ist.
\end{exercise}
752

753
754
755
756
757
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}
 =
758
 x^k_i + \frac{\eta}{A_{ii}} \bigg(b_i-\sum_{j=1}^n A_{ij}x_j^k \bigg)
759
760
761
 \qquad \text{bzw.} \qquad
 x^{k+1} = x^k + \eta D^{-1} ( b - Ax^k).
\end{equation*}
762
763
764
765
766

\subsubsection{Konvergenz}

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

770
771
772
773
774
775
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}.
776

777
Wir greifen die Konvergenztheorie des Jacobi-Verfahrens in Kapitel~\ref{sec:jacobi_smoother_fourier_analysis}
778
wieder auf.
779

780
\subsubsection{Konvergenzgeschwindigkeit}
781
782
783
784
785

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

\bigskip

786
787
Sei $A$ die Matrix des Poisson-Problems auf einem uniformen Gitter
mit Dreiecks- oder Viereckselementen.  Eine Zeile von $A$ ist dann
788
\begin{equation*}
789
    \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.
790
791
\end{equation*}

792
793
794
795
\begin{exercise}
 Rechnen Sie nach dass man tatsächlich diese Matrix erhält.
\end{exercise}

796

797
798
799

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

800
Wir versuchen jetzt, den Spektralradius von $I-CA = I - \eta D^{-1}A$ abzuschätzen:
801
\begin{align*}
802
803
804
805
    \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}}.
806
807
\end{align*}
Der Bruch ist gerade der \emph{Rayleigh-Quotient}
808
809
810
811
812
813
\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.

814
815
\medskip

816
Daraus folgt dass
817
818
819
\begin{align*}
 \rho (I - \eta D^{-1}A )
 & =
820
 \sup \lbrace \abs[\Big]{1- \frac{1}{4} \eta h^2 \lambda} \; : \; \text{$\lambda$ Eigenwert von $A$} \rbrace \\
821
822
 %
 & =
823
 1-2 \eta \sin^2 ( \frac{1}{2} \pi h ).
824
\end{align*}
825
826
Taylor-Entwicklung:
\begin{equation*}
827
    \sin^2 s = s^2 - \frac{1}{3}s^4 + \frac{2}{45}s^6 + \cdots
828
\end{equation*}
829
ergibt
830
\begin{equation*}
831
 \frac{\norm{e^{k+1}}}{\norm{e^k}} \approx \rho (I-D^{-1}A) \approx 1-\frac{\eta}{2} \pi^2h^2
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
\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}
 =
850
 \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,
851
852
853
\end{equation*}
da $\ln x \approx (x-1)-\frac{1}{2} (x-1)^2 + \ldots$ ist.

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

857
858
\medskip

859
Das Jacobi-Verfahren scheint also kein sonderlich gutes Verfahren zu sein.
860
861
862
863
864
865
866

\subsection{Das Gauß-Seidel-Verfahren}

Betrachte noch einmal die Jacobi-Rechenvorschrift:
\begin{equation*}
 x_i^{k+1}
 =
867
868
 \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 ).
869
870
871
\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
872
873
\medskip

874
875
Gauß-Seidel-Verfahren:
\begin{align*}
Sander, Oliver's avatar
Sander, Oliver committed
876
877
 x_i^{k+1}
 & =
878
 \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
879
 & =
880
 \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 ).
881
882
883
\end{align*}
$x_i^{k+1}$ hängt von $x_{i-1}^{k+1}$ ab $\implies$ keine Parallelisierung möglich.

884
885
886
887
888
\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}
889
890
891
892
893
894
895
896
897

\subsubsection{Konvergenz}

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

\medskip

Und in der Tat:
\begin{theorem}
898
899
900
901
Das Gauß-Seidel-Verfahren konvergiert, wenn
\begin{itemize}
 \item $A$ symmetrisch und positiv definit ist und/oder
 \item $A$ irreduzibel und diagonaldominant ist.
902
903
904
905
906
907
\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
908
 \item Es gilt $\rho(I-(D+L)^{-1}A) = ( \rho (I-D^{-1}A ) )^2$.
909
910
911
912
  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}.

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

915
916
917
918
919
  \begin{exercise}
   Finden Sie die genaue Stelle, und sagen Sie mir bescheid.
  \end{exercise}


920
921
922
 %
 \item Deshalb
  \begin{equation*}
Sander, Oliver's avatar
Sander, Oliver committed
923
   \rho(I-(D+L)^{-1}A) = \cos^2(\pi h)
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
  \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
944
  \frac{- \ln R}{\ln (1-\pi^2h^2 )} \approx \frac{\ln R}{\pi^2h^2}
945
946
947
948
949
950
 \end{equation*}
 Iterationen.
\end{itemize}
\end{example}
\emph{Faustregel:} Das Gauß-Seidel-Verfahren braucht nur etwa halb so viele Iterationen wie das Jacobi-Verfahren.

951
952
953
954
955
956
957
958
959
\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}
960
961
962
963
964
965
966
967
968

\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}

969
\subsubsection{Das Verfahren der konjugierten Gradienten (CG-Verfahren)}
970
971

\begin{itemize}
972
973
974
975
976
977
 \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.

978
979
980
981
 \item Geriet zwischenzeitlich in Vergessenheit
 \item Erlebte Revival als iterative Methode
\end{itemize}

982
983
\begin{theorem}[{\cite{shewchuk:1994}}]
Sei $\kappa$ die Kondition der Matrix $A$. Dann gilt nach $k$ Schritten des CG-Verfahrens
984
\begin{align*}
985
 \norm{e^k}_A \leq 2 ( \frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1} )^k \cdot \norm{e^0}_A.
986
\end{align*}
987
988
\end{theorem}

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

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

994
995
 \item Aber man benötigt wieder einen Vorkonditionierer, sonst ist auch
   das CG-Verfahren zu langsam.
996

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

999
1000
 \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.
1001

1002
1003
 \item Die Konstruktion passender $C$ bildet ein großes Gebiet innerhalb
   der Numerik.
1004

1005
 \item Eine Möglichkeit: Mehrgitter!
1006
1007
\end{itemize}

1008
\section{Mehrgitterverfahren: Der anschauliche Zugang}
1009

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

1013
\subsection{Glättung}
1014

1015
Gauß--Seidel- und Jacobi-Verfahren sind für typische Finite-Elemente-Matrizen
1016
extrem langsam.
1017

1018
1019
1020
1021
\medskip

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

1023
1024
1025
\smallskip

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

1027
1028
1029
\bigskip

\emph{Beispiel:} Wir betrachten wieder das Poisson-Problem
1030
1031
1032
1033
1034
\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$.
1035
1036
1037

\medskip

1038
1039
1040
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:
1041
1042
1043
1044
1045
1046
1047
\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.

1048
\medskip
1049
1050
1051

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

1052
\medskip
1053
1054
1055
1056

Hier sind einige der ersten Iterierten:

\begin{center}
1057
 \begin{tabular}{ccc}
1058
1059

  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_zero_iterate0} &
1060
1061
  \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} \\
1062

1063
  $x^0$ & $x^1$ & $x^2$ \\
1064

1065
  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_zero_iterate3} &
1066
1067
1068
  \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} \\

1069
  $x^3$ & $x^{10}$ & $x^{50}$
1070
1071
1072
1073
1074
1075
1076
1077
 \end{tabular}
\end{center}

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

\bigskip

1078
Jetzt zeigen wir was wir mit \glqq Gauß--Seidel glättet\grqq{} meinen.
1079
1080
1081

\medskip

1082
1083
1084
1085
Wir lösen nochmal die Gleichung, fangen jetzt aber von einem
verrauschten Startwert an.

\begin{center}
1086
 \begin{tabular}{ccc}
1087
  \includegraphics[width=0.3\textwidth,trim=0 50 0 0]{gs_from_noise_iterate0} &
1088
1089
  \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} \\
1090

1091
  $x^0$ & $x^1$ & $x^2$ \\
1092

1093
  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_noise_iterate3} &
1094
1095
1096
  \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} \\

1097
  $x^3$ & $x^{10}$ & $x^{50}$
1098
1099
1100
 \end{tabular}
\end{center}

1101
\begin{itemize}
1102
 \item Das Verfahren konvergiert immer noch.
1103

1104
 \item Es scheint weder schneller noch langsamer geworden zu sein.
1105

1106
1107
1108
 \item Das Rauschen ist aber schon nach sehr wenigen Iterationen weitestgehend
   verschwunden.
\end{itemize}
1109

1110
1111
1112
1113
1114
1115
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
1116
   \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).
1117
1118
1119
1120
\end{equation*}
Dabei tauchen in den Summen nur die Nachbarknoten auf.
\begin{itemize}
 \item Zumindest in diesem Modellproblem sind die entsprechenden Koeffizienten
1121
   (fast) alle gleich.
1122
1123
1124
1125
1126
1127

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

\bigskip

1128
So richtig überzeugend ist diese Argumentation trotzdem noch nicht:
1129

1130
1131
1132
\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?
1133

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

1137
1138
In Wirklichkeit ist es etwas anders:

1139
\begin{itemize}
1140
1141
1142
 \item Das Gauß--Seidel-Verfahren glättet \emph{nicht} die Iterierten.

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

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

1147
\subsection{Fourier-Analyse}
1148
\label{sec:jacobi_smoother_fourier_analysis}
1149

1150
1151
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.
1152

1153
1154
\medskip