skript-mehrgitter-sander.tex 37.5 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
14

\usepackage[utf8]{inputenc}

\usepackage{longtable,tabularx,tabulary,booktabs}
\usepackage{array}
15
\usepackage[linesnumbered,inoutnumbered]{algorithm2e}
16
17
\usepackage{caption}
\usepackage{url}
Sander, Oliver's avatar
Sander, Oliver committed
18
\usepackage{overpic}
19
20
21
22
23
24
\usepackage{graphicx}
\usepackage{tikz}

\usepackage{amsfonts}
\usepackage{lmodern}
\usepackage[ngerman]{babel}
25
\usepackage{csquotes}
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{mathtools}
\usepackage{colonequals}
\usepackage{enumitem}

\usepackage[defernumbers=true,
            maxbibnames=99,  % Never write 'et al.' in a bibliography
            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}
\usepackage{todonotes}
\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{\T}{\mathcal T}
\newcommand{\curl}{\operatorname{curl}}
\renewcommand{\div}{\operatorname{div}}
54
\newcommand{\grad}{\operatorname{grad}}
55
56
\newcommand{\tr}{\operatorname{tr}}
\DeclareMathOperator*{\argmin}{arg\,min}
Sander, Oliver's avatar
Sander, Oliver committed
57
\newcommand{\Gammatight}[1]{{\Gamma\hspace{-0.8mm}_{#1}}}  % A \Gamma with a subscript, and extra kerning
58
59
60
61
62

% Bold letters
\newcommand{\bn}{\mathbf{n}}
\newcommand{\bv}{\mathbf{v}}

63
64
\DeclarePairedDelimiter{\abs}{\lvert}{\rvert}
\DeclarePairedDelimiter{\norm}{\lVert}{\rVert}
65
66
67
68
69
70
71
72
73
74
75
76
77

\theoremstyle{plain} %Text ist Kursiv
\newtheorem{theorem}{Satz}[chapter]
\newtheorem{hilfssatz}{Hilfssatz}[chapter]
\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}
78
\newtheorem{exercise}[theorem]{Übung}
79
80
81
82
83
84
85
86
87
88
89
90
91
92

\swapnumbers

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

\begin{document}

\begin{titlepage}

\begin{center}

\vspace*{0.3\textheight}

93
{\Huge Mehrgitterverfahren}
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
\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
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
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
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
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
\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$.
Wir betrachten als Beispiel die {\em Poisson-Gleichung}
\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}
 a(u,v) = l(v)
\end{equation}
mit
\begin{equation*}
 a(v,w) = \int_\Omega \nabla v \nabla w \, dx
\quad \text{und} \quad
 l(v) = \int_{\Gamma_N} g v \, ds
+
\int_\Omega f v \, dx.
\end{equation*}
Man nennt \eqref{eq:schwache_formulierung} die {\em schwache Formulierung}
(oder {\em Variationsformulierung})
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
a(u,v) = l(v)
\qquad v \in H
\end{equation*}
for jedes $l \in H'$ eine eindeutig bestimmte Lösung.
\end{theorem}




\section{Diskretisierung (Galerkin-Verfahren)}

Sei $V_h$ ein endlichdimensionaler Teilraum von $H^1_{\Gammatight{D}}$ (mit Dimension $n$).
Der Index $h$
bezeichne einen Ordnungsparameter.  Die {\em diskrete Variationsformulierung}
lautet
\begin{equation}
\label{eq:diskrete_formulierung}
 u_h \in V_h \quad : \quad a(u_h,v_h) = l(v_h)
\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}
 u_h \in V_h \quad : \quad a(u_h, \varphi_i) = l(\varphi_i)
\qquad \text{für alle $0 \le i < n$}.
\end{equation*}
Dies wiederum entspricht dem linearen Gleichungssystem
267
268
\begin{equation}
\label{eq:fe_lineares_gleichungssystem}
Sander, Oliver's avatar
Sander, Oliver committed
269
 A\bar{u} = b,
270
\end{equation}
Sander, Oliver's avatar
Sander, Oliver committed
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
wobei
\begin{alignat*}{2}
A &\in \R^{n \times n},
\qquad &A_{ij} & = a(\varphi_i, \varphi_j)
=
\int_\Omega \nabla \varphi_i \nabla \varphi_j \, dx \\
%
b &\in \R^n,
\qquad &
b_i & = l(\varphi_i)
=
\int_{\Gamma_N} g \varphi_i \, ds
+
\int_\Omega f \varphi_i \, dx
\end{alignat*}
und $\bar{u} \in \R^n$ die Koeffizienten von $u_h$ bzgl.\ der Basis $\{\varphi_i\}$ sind.

Die Matrix $A$ nennt man {\em Steifigkeitsmatrix}.  Dieser Ausdruck kommt aus den
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$
nennt man deshalb manchmal auch {\em Lastvektor}.

\emph{Beachte:} Ohne weitere Annahmen ist an dieser Stelle davon auszugehen, dass $A$
voll besetzt ist.  Das bedeutet, dass fast jeder Eintrag von $A$ von Null verschieden ist,
wobei der Ausdruck \glqq fast jeder\grqq{} hier nicht mathematisch streng gemeint ist.

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

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

Sei $\Omega$ von jetzt ab {\em polygonal berandet}.  Wir füllen
$\Omega$ mit einer Triangulierung / einem {\em Gitter}.  Ein Beispielgitter für
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}$
von abgeschlossenen Dreiecken $t$ heißt {\em Triangulierung} von $\Omega$, falls gilt
\begin{enumerate}
 \item Die Menge der Dreiecke überdeckt den Abschluss des Gebiets
\begin{equation*}
 \overline{\Omega} = \bigcup_{t \in \mathcal{T}} t.
\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*}
 h = \max_{t \in \mathcal{T}} \operatorname{diam} t,
\end{equation*}
der als Index in $V_h$ auftaucht, beschreibt traditionell die Feinheit des Gitters.
Oft nennt man $h$ auch {\em Gitterweite}.

\begin{definition}
Der Raum
\begin{equation*}
 V_h^{(1)}
\colonequals
\{ v \in C(\overline{\Omega})
\; | \;
\text{$v$ ist affin auf $t$ für alle $t \in \mathcal{T}$, und $v|_{\Gamma_D} = 0$} \}
\end{equation*}
heißt Raum der Lagrange-Elemente erster Ordnung
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}
\caption{Links: eine Finite-Elemente Funktion.  Rechts: ein Element der Knotenbasis}
\label{fig:beispiel-fe-funktion}
\end{figure}

Als Basis von $V_h^{(1)}$ wählt man gewöhnlich die {\em Knotenbasis}.  Sei $\mathcal V$
die Menge der Knoten des Gitters.  Dann ist die Knotenbasis die Menge
der Funktionen $\varphi_i \in V_h^{(1)}$ mit
\begin{equation*}
 \varphi_i (v_j) = \delta_{ij}
\qquad
\text{für alle $v_j \in \mathcal{V}$}.
\end{equation*}
Es ist also insbesondere $\dim V_h^{(1)} = \abs{\mathcal V \setminus \mathcal{V}_{\Gamma_D}}$
(Abbildung~\ref{fig:beispiel-fe-funktion}).

{\em Wichtig:} Da die Basisfunktionen $\varphi_i$ einen `lokalen'
Träger haben, so sind die meisten Matrixeinträge
\begin{equation*}
 A_{ij} = a(\varphi_i, \varphi_j)
=
\int_\Omega \nabla \varphi_i \nabla \varphi_j \, dx
\end{equation*}
gleich Null!  Die Matrix $A$ ist also {\em dünnbesetzt}.
Das ist eine wichtige Eigenschaft, denn müsste man alle Einträge von $A$ in Betracht ziehen,
so wären das \glqq zu viele\grqq{} (nämlich quadratisch viele).
Nur bei einer dünnbesetzten Matrix bleibt der Aufwand handhabbar.
Denn $n$ kann sehr groß werden; aktuell sind Werte im Bereich weniger Millionen auf einem
einzelnen Rechner realistisch.  Auf einem großen Parallelrechner können es durchaus einige
Milliarden werden.

\begin{theorem}
 Es sei $h$ klein genug, $u \in H^1_0(\Omega) \cap H^{2}(\Omega)$,
$\Omega \in \R^2$.  Dann gilt die a priori Abschätzung
\begin{equation*}
 \norm{u - u_h}_1 \le C h \abs{u}_{2}
\end{equation*}
mit einer Konstanten $C$.
\end{theorem}

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

399
\chapter{Mehrgitterverfahren}
Sander, Oliver's avatar
Sander, Oliver committed
400

401
\section{Iterative Lösungsverfahren für dünnbesetzte lineare Gleichungssysteme}
402

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

406
\subsection{Eigenschaften von Finite-Elemente-Matrizen}
407
408
409
\subsubsection{Größe}
Die Matrizen aus dem obigen Beispiel können sehr groß werden.
\begin{itemize}
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
 \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.

 \item Eine Zahl in doppelter Genauigkeit braucht $8$~Byte.  In ein GB RAM
   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}

Die Matrizen sind schlecht konditioniert.

\emph{Beispiel:} $\Omega = [0,1]^2$, strukturiertes Dreiecksgitter mit Kantenlänge $h$.
\begin{itemize}
 \item Die Eigenwerte von $A$ sind~\cite[Kapitel~12.3.3]{dahmen_reusken:2008}
  \todo[inline]{Kommt noch}

 \item Die Kondition ist deshalb (die Matrix ist symmetrisch)
  \begin{equation*}
   \kappa
   =
   \frac{\lambda_\text{max}}{\lambda_\text{min}}
   =
  \end{equation*}
  \todo[inline]{Kommt noch}

 \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}
454
455
456
\end{itemize}


457
458
459
460
461

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

462
463

\subsubsection{Dünnbesetztheit}
464
465
466

Mehr zu lesen findet man z.B. bei \citet{pissanetzky:1984}.

467
\begin{itemize}
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
\item Ein Matrixeintrag $A_{ij} = a(\varphi_i, \varphi_j)$ ist genau dann
  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

\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.
484
485
486
487
488
489
490
491
492
493
494
495
496
497
\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$\\
498
     \For{alle Spalten $j$ in denen $A_{ij} \neq 0$}
499
500
501
502
503
     {
       $v_i = v_i + A_{ij} w_j$
     }
   }
  \end{algorithm}
504
Das braucht nur $O(\# \text{Nichtnulleinträge})$ Operationen (statt $O(n^2)$).
505

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

508
\emph{Erinnerung:} Gauß-Elimination braucht $O(n^3)$ Rechenoperationen!
509

510
\subsection{Direkte Verfahren für dünnbesetzte Gleichungen}
511

512
513
\begin{itemize}
 \item Es gibt direkte Verfahren zum Lösen von dünnbesetzten Gleichungssystemen.
514

515
 \item Sehr kompliziert / sehr ausgefeilt.
516

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

519
520
521
522
523
524
525
526
527
528
529
530
 \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}
531
532
533
534


\subsection{Lineare iterative Verfahren}

535
536
537
538
539
540
Iterative Verfahren konstruieren eine Folge $\{ x^k \}_{k \in \N}$,
die (hoffentlich) gegen die Lösung $x^*$ von $Ax=b$ konvergiert.

\medskip

Bei linearen iterativen Verfahren ist die Abbildung
541
\begin{equation*}
542
543
544
 e^k \colonequals x^k - x^*
 \qquad \mapsto \qquad
 e^{k+1} \colonequals x^{k+1} - x^*
545
\end{equation*}
546
linear.
547
548
549

\bigskip

550
Der normale Ansatz ist~\cite{dahmen_reusken:2008}:
551
552
553
554
555

Schreibe $Ax=b$ als Fixpunktgleichung
\begin{equation*}
	x=x+C(b-Ax)
\end{equation*}
556
557
mit $C \in \R^{n \times n}$ nichtsingulär. Die Größe $r(x) \colonequals b - Ax$
nennt man \emph{Residuum}.  Die Matrix $C$ heißt \emph{Vorkonditionierer}.
558
559
560

Dafür machen wir jetzt eine Fixpunktiteration:
\begin{equation*}
561
  x^{k+1}\colonequals x^k+C (b-Ax^k)=(I-CA)x^k+Cb,
562
563
564
565
566
567
568
569
570
571
572
573
  \qquad
  k = 0, 1, 2, \dots
\end{equation*}

Unter welchen Umständen konvergiert dieses Verfahren?

\medskip

Der Fehler im $k$-ten Iterationsschritt ist $e^k=x^k-x^*$. Es gilt
\begin{equation*}
 e^{k+1}=x^{k+1}-x^*
 =
574
\underbrace{(I-CA)}_{\text{Iterationsmatrix}}e^k.
575
576
577
578
579
580
581
582
583
584
585
586
587
588
\end{equation*}
Also gilt
\begin{equation}
\label{equa:eqerrlin}
 e^k=(I-CA)^ke^0
 \qquad
 \forall k=0,1,2,\ldots
\end{equation}

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


589
Die Fehlerfortpflanzung \eqref{equa:eqerrlin} ist tatsächlich linear.
590
591
592
593

\bigskip

Für die Konvergenz gilt der folgende wichtige Satz.
594
595

Für jede Vektornorm mit assoziierter Matrixnorm erhält man
596
\begin{theorem}
597
598
 Das Verfahren konvergiert für jeden Startwert $x^0 \in \R^n$ gegen die Lösung
 von $Ax = b$, wenn $\norm{I - CA} < 1$.
599
600
\end{theorem}

601
Da $\rho (B) \leq \norm{B}$ für jede submultiplikative Matrixnorm gilt sogar
602
603
604
605
606
607
608
609
610
611
612
613
614

[Denn: Sei $v$ ein Eigenvektor von $B$ zum Eigenwert $\lambda$.  Dann ist
\begin{equation*}
 \abs{\lambda}\norm{v}
 =
 \norm{\lambda v}
 =
 \norm{B v}
 \le
 \norm{B}\norm{v}.
\end{equation*}
Deshalb gilt $\abs{\lambda} \le \norm{B}$ für alle Eigenwerte $\lambda$ von $B$.]

615
616
617
\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}
618

619
Je nach Wahl des Vorkonditionierers $C$ erhält man unterschiedliche Verfahren:
620
\begin{itemize}
621
 \item Das Jacobi-Verfahren
622
  \begin{equation*}
623
624
   C = D^{-1}
  \end{equation*}
625

626
627
628
629
 \item Das Gauß--Seidel-Verfahren
  \begin{equation*}
   C = (D+L)^{-1}
  \end{equation*}
630

631
 \item Das symmetrische Gauß--Seidel-Verfahren
632
  \begin{equation*}
633
   C = (D-L)^{-1}+(D-R)^{-1}-(D-R)^{-1}A(D-L)^{-1}
634
635
636
  \end{equation*}
\end{itemize}

637
Es ergibt sich folgendes Dilemma:
638
639
640
641
642
643
644
645
\begin{enumerate}
    \item $C$ soll $A^{-1}$ möglichst gut approximieren
    \item Die Operation $y \mapsto Cy$ soll möglichst billig sein
\end{enumerate}
\emph{Beachte:} Die Matrix $C$ wird nie explizit ausgerechnet!


\subsection{Das Jacobi-Verfahren}
646
647
648
649

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

650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
Wir nehmen im Folgenden an, dass $a_{ii} \neq 0$ für alle $i=1,\ldots,n$.

\medskip

Betrachte das lineare Gleichungssystem:
\begin{align*}
    a_{11} x_1+a_{12} x_2+\ldots +a_{1m} x_m & = b_1 \\
    a_{12} x_1+a_{22} x_2+\ldots +a_{2m} x_m & = b_2 \\
    & \vdots
\end{align*}
\emph{Idee:} Löse $i$-te Zeile nach $x_i$ auf für $i=1,\ldots,n$.
\begin{align*}
    x_1 & =\frac{1}{a_{11}} \left(b_1-a_{12} x_2-\ldots -a_{1n} x_n \right) \\
    x_2 & =\frac{1}{a_{22}} \left(b_2-a_{21} x_1-\ldots -a_{2n} x_n \right) \\
    & \vdots
\end{align*}
Mache daraus ein iteratives Verfahren:
\begin{align*}
    x_1^{k+1} & =\frac{1}{a_{11}} \left(b_1-a_{12}x_2^k-\ldots -a_{1n}x_n^k \right) \\
    x_2^{k+1} & =\frac{1}{a_{22}} \left(b_2-a_{21}x_1^k-\ldots -a_{2n}x_n^k \right) \\
    & \vdots
\end{align*}
Für $i=1,\ldots,n$
\begin{equation*}
674
675
676
677
678
 x_i^{k+1}
 =
 \frac{1}{a_{ii}} \bigg(b_i-\sum_{\substack{j=1\\j\neq i}}^n a_{ij}x_i^k \bigg)
 =
 x^k_i + \frac{1}{a_{ii}} \bigg(b_i-\sum_{j=1}^n a_{ij}x_i^k \bigg).
679
\end{equation*}
680
681
682
Beachte:
\begin{itemize}
 \item Die Rechnungen für die verschiedenen $i$ sind voneinander unabhängig $\implies$ leicht zu parallelisieren.
683

684
685
 \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}
686

687
688
689
690
\begin{exercise}
 Rechnen Sie nach dass dies wirklich ein lineares Verfahren mit Vorkonditionierer
 $C = D^{-1}$ ist.
\end{exercise}
691

692
\todo[inline]{Hier muss direkt das gedämpfte Verfahren eingeführt werden.}
693
694
695
696
697
698
699
700
701

\subsubsection{Konvergenz}

Wir wenden das bekannte Konvergenzkriterium an:
\begin{theorem}
Das Jacobi Verfahren konvergiert genau dann, wenn $\rho \left(I-D^{-1}A \right) <1$.
\end{theorem}
Leider gilt diese Bedingung nicht immer.

702
Ein paar direkte Charakterisierungen findet man z.B.\ bei \citet{dahmen_reusken:2008}.
703

704
705
Wir greifen die Konvergenztheorie des Jacobi-Verfahrens in Kapitel~\ref{sec:jacobi_smoother}
wieder auf.
706

707
\subsubsection{Konvergenzgeschwindigkeit}
708
709
710
711
712

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

\bigskip

713
714
Sei $A$ die Matrix des Poisson-Problems auf einem uniformen Gitter
mit Dreiecks- oder Viereckselementen.  Eine Zeile von $A$ ist dann
715
716
717
718
\begin{equation*}
    \frac{1}{h^2} \left(4u_{i,j}-u_{i-1,j}-u_{i+1,j}-u_{i,j-1}-u_{i,j+1} \right)=f_i.
\end{equation*}

719
720
721
722
\begin{exercise}
 Rechnen Sie nach dass man tatsächlich diese Matrix erhält.
\end{exercise}

723

724
725
726
727
728
729
730
731
732
733
734

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

Wir versuchen jetzt, den Spektralradius von $I-CA = I - D^{-1}A$ abzuschätzen:
\begin{align*}
    \rho (I-D^{-1}A)
    & = \sup_{x \neq 0} \abs[\bigg]{\frac{x^T \left(I-D^{-1}A \right) x}{\norm{x}^2}} \\
    & = \sup_{x \neq 0} \abs[\bigg]{1-\frac{x^T D^{-1}Ax}{\norm{x}^2}} \\
    & = \sup_{x \neq 0} \abs[\bigg]{1-\frac{1}{4}h^2 \frac{x^TAx}{\norm{x}^2}}.
\end{align*}
Der Bruch ist gerade der \emph{Rayleigh-Quotient}
735
736
737
738
739
740
741
742
\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.

Daraus folgt dass
\begin{equation*}
743
744
745
 \rho (I-D^{-1}A )
 =
 \sup \left\lbrace \abs[\Big]{1- \frac{1}{4} h^2 \lambda} \; : \; \text{$\lambda$ Eigenwert von $A$} \right\rbrace.
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
\end{equation*}

Für die spezielle Matrix können wir die Eigenwerte ausrechnen.  Diese sind
alle nichtnegativ.

Der größte Eigenwert von $A$ ist~\cite[Kapitel~12.3.3]{dahmen_reusken:2008}:

\begin{equation*}
 \lambda =\frac{8}{h^2} \sin^2 \left( \frac{1}{2} \pi h \right)
\end{equation*}
Es folgt, dass
\begin{equation*}
 \rho (I-D^{-1}A)=1-2 \sin^2 \left( \frac{1}{2} \pi h \right)=\cos (\pi h)
\end{equation*}
Taylor-Entwicklung:
\begin{equation*}
    \cos \left( \pi h \right)=1-\frac{1}{2} \pi^2 h^2+\ldots
\end{equation*}
Also ist
\begin{equation*}
 \frac{\norm{e^{k+1}}}{\norm{e^k}} \approx \rho (I-D^{-1}A) \approx 1-\frac{1}{2} \pi^2h^2
\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}
 =
 \frac{-\ln R}{\ln \rho \left(I-D^{-1}A \right)} \approx \frac{-\ln R}{\ln \left(1-\frac{1}{2}\pi^2 h^2 \right)} \approx \frac{2}{\pi^2 h^2} \ln R,
\end{equation*}
da $\ln x \approx (x-1)-\frac{1}{2} (x-1)^2 + \ldots$ ist.

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

792
Das Jacobi-Verfahren scheint also kein sonderlich gutes Verfahren zu sein.
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811

\subsection{Das Gauß-Seidel-Verfahren}

Betrachte noch einmal die Jacobi-Rechenvorschrift:
\begin{equation*}
 x_i^{k+1}
 =
 \frac{1}{a_{ii}} \left(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 \right).
\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}$.

Gauß-Seidel-Verfahren:
\begin{align*}
    x_i^{k+1} & =\frac{1}{a_{ii}} \left(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 \right) \\
    & = \frac{1}{a_{ii}} \left(b_i-\sum_{j=1}^{i-1} a_{ij}x_j^{k+1} -\sum_{j=i+1}^n a_{ij} x_j^k \right).
\end{align*}
$x_i^{k+1}$ hängt von $x_{i-1}^{k+1}$ ab $\implies$ keine Parallelisierung möglich.

812
813
814
815
816
\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}
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839

\subsubsection{Konvergenz}

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

\medskip

Und in der Tat:
\begin{theorem}
Das Gauß-Seidel-Verfahren konvergiert, wenn \begin{itemize}
	\item $A$ symmetrisch und positiv definit ist und/oder
	\item $A$ irreduzibel und diagonal dominant ist.
\end{itemize}
\end{theorem}
\subsubsection{Konvergenzgeschwindigkeit}
\begin{example}
Sei $A$ die Matrix des Poisson-Problems für ein quadratisches Gebiet.
\begin{itemize}
 \item Es gilt $\rho(I-(D-L)^{-1}A) = ( \rho (I-D^{-1}A ) )^2$.
  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}.

840
841
842
843
844
  \begin{exercise}
   Finden Sie die genaue Stelle, und sagen Sie mir bescheid.
  \end{exercise}


845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
 %
 \item Deshalb
  \begin{equation*}
   \rho(I-(D-L)^{-1}A) = \cos^2(\pi h)
  \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
  \frac{- \ln R}{\ln \left(1-\pi^2h^2 \right)} \approx \frac{\ln R}{\pi^2h^2}
 \end{equation*}
 Iterationen.
\end{itemize}
\end{example}
\emph{Faustregel:} Das Gauß-Seidel-Verfahren braucht nur etwa halb so viele Iterationen wie das Jacobi-Verfahren.

876
877
878
879
880
881
882
883
884
\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}
885
886
887
888
889
890
891
892
893

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

894
\subsubsection{Das Verfahren der konjugierten Gradienten (CG-Verfahren)}
895
896

\begin{itemize}
897
898
899
900
901
902
 \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.

903
904
905
906
 \item Geriet zwischenzeitlich in Vergessenheit
 \item Erlebte Revival als iterative Methode
\end{itemize}

907
908
\begin{theorem}[{\cite{shewchuk:1994}}]
Sei $\kappa$ die Kondition der Matrix $A$. Dann gilt nach $k$ Schritten des CG-Verfahrens
909
\begin{align*}
910
 \norm{e^k}_A \leq 2 \left( \frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1} \right)^k \cdot \norm{e^0}_A.
911
\end{align*}
912
913
\end{theorem}

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

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

919
920
 \item Aber man benötigt wieder einen Vorkonditionierer, sonst ist auch
   das CG-Verfahren zu langsam.
921

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

924
925
 \item $C$ muss so gewählt werden dass $\kappa(CA)$ klein ist, aber $Cv$, $v \in \R^n$
   billig ausgerechnet werden kann.
926

927
928
 \item Die Konstruktion passender $C$ bildet ein großes Gebiet innerhalb
   der Numerik.
929

930
 \item Eine Möglichkeit: Mehrgitter!
931
932
\end{itemize}

933
\section{Mehrgitterverfahren: Der anschauliche Zugang}
934

935
\subsection{Glättung}
936

937
938
Gauß--Seidel-Verfahren und Jacobi-Verfahren sind für typische Finite-Elemente-Matrizen
extrem langsam.
939

940
941
942
943
944
945
946
\medskip

Sie haben allerdings eine wichtige Eigenschaft, die wir uns
für Mehrgitterverfahren zunutze machen:
\begin{itemize}
 \item Die zwei Verfahren \emph{glätten}.
\end{itemize}
947

948
\emph{Beispiel:} Wir betrachten das Poisson-Problem
949

950
951
952
953
954
\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$.
955
956
957

\medskip

958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
Hier ist die Finite-Elemente-Lösung auf einem strukturierten Dreiecksgitter
mit $16 \times 16 \times 2$ Elementen, dargestellt als Höhenplot:
\begin{center}
 \includegraphics[width=0.6\textwidth,trim=0 50 0 10]{gs_from_zero_iterate50}
\end{center}

Wir versuchen jetzt, das dazugehörige algebraische Problem $Ax = b$
mit dem Gauß--Seidel-Verfahren zu lösen.

\bigskip

Eine Iteration ist:

\begin{algorithm}[H]
 \SetAlgoLined
  % \caption{...}
 \For{alle Zeilen $i$}
 {
  \begin{equation*}
   x^{k+1}_i
    \longleftarrow
   \frac{1}{A_{ii}}
     \bigg( b_i - \sum_{j=1}^{i-1} A_{ij} x^{k+1}_j - \sum_{j=1+1}^{n} A_{ij} x^k_j \bigg)
  \end{equation*}
 }
\end{algorithm}

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

\bigskip

Hier sind einige der ersten Iterierten:

\begin{center}
 \begin{tabular}{cc}

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

  iterate $x^0$ & iterate $x^1$ \\

  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_zero_iterate2} &
  \includegraphics[width=0.3\textwidth,trim=0 100 0 0]{gs_from_zero_iterate3} \\
For faster browsing, not all history is shown. View entire blame