skript-mehrgitter-sander.tex 44.8 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
\usepackage{graphicx}
\usepackage{tikz}
21
\usepackage{pgfplots}
22
23
24
25

\usepackage{amsfonts}
\usepackage{lmodern}
\usepackage[ngerman]{babel}
26
\usepackage{csquotes}
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
54
\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}}
55
\newcommand{\grad}{\operatorname{grad}}
56
57
\newcommand{\tr}{\operatorname{tr}}
\DeclareMathOperator*{\argmin}{arg\,min}
Sander, Oliver's avatar
Sander, Oliver committed
58
\newcommand{\Gammatight}[1]{{\Gamma\hspace{-0.8mm}_{#1}}}  % A \Gamma with a subscript, and extra kerning
59
60
61
62
63

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

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

\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}
79
\newtheorem{exercise}[theorem]{Übung}
80
81
82
83
84
85
86
87
88
89
90
91
92
93

\swapnumbers

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

\begin{document}

\begin{titlepage}

\begin{center}

\vspace*{0.3\textheight}

94
{\Huge Mehrgitterverfahren}
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
129
\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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
\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$.
153
Wir betrachten als Beispiel die \emph{Poisson-Gleichung}
Sander, Oliver's avatar
Sander, Oliver committed
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
\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}
204
 a(u,v) = \ell(v)
Sander, Oliver's avatar
Sander, Oliver committed
205
206
207
\end{equation}
mit
\begin{equation*}
208
209
210
 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
211
212
213
+
\int_\Omega f v \, dx.
\end{equation*}
214
215
Man nennt \eqref{eq:schwache_formulierung} die \emph{schwache Formulierung}
(oder \emph{Variationsformulierung})
Sander, Oliver's avatar
Sander, Oliver committed
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
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
237
a(u,v) = \ell(v)
Sander, Oliver's avatar
Sander, Oliver committed
238
239
\qquad v \in H
\end{equation*}
240
for jedes $\ell \in H'$ eine eindeutig bestimmte Lösung.
Sander, Oliver's avatar
Sander, Oliver committed
241
242
243
244
245
246
247
248
249
\end{theorem}




\section{Diskretisierung (Galerkin-Verfahren)}

Sei $V_h$ ein endlichdimensionaler Teilraum von $H^1_{\Gammatight{D}}$ (mit Dimension $n$).
Der Index $h$
250
bezeichne einen Ordnungsparameter.  Die \emph{diskrete Variationsformulierung}
Sander, Oliver's avatar
Sander, Oliver committed
251
252
253
lautet
\begin{equation}
\label{eq:diskrete_formulierung}
254
 u_h \in V_h \quad : \quad a(u_h,v_h) = \ell(v_h)
Sander, Oliver's avatar
Sander, Oliver committed
255
256
257
258
259
260
261
262
263
\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}
264
 u_h \in V_h \quad : \quad a(u_h, \varphi_i) = \ell(\varphi_i)
Sander, Oliver's avatar
Sander, Oliver committed
265
266
267
\qquad \text{für alle $0 \le i < n$}.
\end{equation*}
Dies wiederum entspricht dem linearen Gleichungssystem
268
269
\begin{equation}
\label{eq:fe_lineares_gleichungssystem}
270
 Ax = b,
271
\end{equation}
Sander, Oliver's avatar
Sander, Oliver committed
272
273
274
wobei
\begin{alignat*}{2}
A &\in \R^{n \times n},
275
\qquad &A_{ij} & \colonequals a(\varphi_i, \varphi_j)
Sander, Oliver's avatar
Sander, Oliver committed
276
277
278
279
280
=
\int_\Omega \nabla \varphi_i \nabla \varphi_j \, dx \\
%
b &\in \R^n,
\qquad &
281
b_i & \colonequals \ell(\varphi_i)
Sander, Oliver's avatar
Sander, Oliver committed
282
283
284
285
286
=
\int_{\Gamma_N} g \varphi_i \, ds
+
\int_\Omega f \varphi_i \, dx
\end{alignat*}
287
und $x \in \R^n$ die Koeffizienten von $u_h$ bzgl.\ der Basis $\{\varphi_i\}$ sind.
Sander, Oliver's avatar
Sander, Oliver committed
288

289
290
291
\bigskip

Die Matrix $A$ nennt man \emph{Steifigkeitsmatrix}.  Dieser Ausdruck kommt aus den
Sander, Oliver's avatar
Sander, Oliver committed
292
293
294
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$
295
nennt man deshalb manchmal auch \emph{Lastvektor}.
Sander, Oliver's avatar
Sander, Oliver committed
296
297
298
299
300
301
302
303


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

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

304
305
306
307
\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
308
309
310
311
312
313
314
315
316
317
318
319
320
321
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}$
322
von abgeschlossenen Dreiecken $T$ heißt \emph{Triangulierung} von $\Omega$, falls gilt
Sander, Oliver's avatar
Sander, Oliver committed
323
324
325
\begin{enumerate}
 \item Die Menge der Dreiecke überdeckt den Abschluss des Gebiets
\begin{equation*}
326
 \overline{\Omega} = \bigcup_{T \in \mathcal{T}} T.
Sander, Oliver's avatar
Sander, Oliver committed
327
328
329
330
331
332
333
334
\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*}
335
 h \colonequals \max_{T \in \mathcal{T}} \operatorname{diam} T,
Sander, Oliver's avatar
Sander, Oliver committed
336
337
\end{equation*}
der als Index in $V_h$ auftaucht, beschreibt traditionell die Feinheit des Gitters.
338
Oft nennt man $h$ auch \emph{Gitterweite}.
Sander, Oliver's avatar
Sander, Oliver committed
339
340
341
342

\begin{definition}
Der Raum
\begin{equation*}
343
 V_h^{p}
Sander, Oliver's avatar
Sander, Oliver committed
344
345
346
\colonequals
\{ v \in C(\overline{\Omega})
\; | \;
347
\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
348
\end{equation*}
349
heißt Raum der Lagrange-Elemente $p$-Ordnung
Sander, Oliver's avatar
Sander, Oliver committed
350
351
352
353
354
355
356
357
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}
358
\caption{Links: eine Finite-Elemente Funktion.  Rechts: ein Element der Knotenbasis erster Ordnung}
Sander, Oliver's avatar
Sander, Oliver committed
359
360
361
\label{fig:beispiel-fe-funktion}
\end{figure}

362
363
364
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
365
366
367
\begin{equation*}
 \varphi_i (v_j) = \delta_{ij}
\qquad
368
\text{für alle $v_j \in \mathcal{L}$}.
Sander, Oliver's avatar
Sander, Oliver committed
369
\end{equation*}
370
Es ist also insbesondere $\dim V_h^{p} = \abs{\mathcal L \setminus \mathcal{L}_{\Gamma_D}}$
Sander, Oliver's avatar
Sander, Oliver committed
371
372
(Abbildung~\ref{fig:beispiel-fe-funktion}).

373
374
375
376
\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
377
378

\begin{theorem}
379
380
381
 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
382
\begin{equation*}
383
 \norm{u - u_h}_1 \le C h^p \abs{u}_{p+1}
Sander, Oliver's avatar
Sander, Oliver committed
384
385
386
387
388
389
390
391
\end{equation*}
mit einer Konstanten $C$.
\end{theorem}

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

392
\chapter{Mehrgitterverfahren}
Sander, Oliver's avatar
Sander, Oliver committed
393

394
\section{Iterative Lösungsverfahren für dünnbesetzte lineare Gleichungssysteme}
395

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

399
\subsection{Eigenschaften von Finite-Elemente-Matrizen}
400
401
402

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

403
\subsubsection{Größe}
404
Finite-Elemente-Matrizen können sehr groß werden:
405
\begin{itemize}
406
407
408
409
410
411
412
 \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.

413
 \item Eine Zahl in doppelter Genauigkeit braucht $8$~Byte.  In ein Gigabyte Hauptspeicher
414
415
416
417
418
419
420
421
422
   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}

423
424
Die Matrizen zu elliptischen Problemen sind positiv definit, aber sie
sind schlecht konditioniert.
425

426
\medskip
427

428
429
430
431
432
433
434
435
436
437
438
439
440
\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)
441
  \begin{equation*}
442
   \kappa(A)
443
444
445
   =
   \frac{\lambda_\text{max}}{\lambda_\text{min}}
   =
446
447
448
449
450
451
452
   \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).
453
454
455
456
457
458
459
460
461
462
463
464
465
  \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}
466
467
\end{itemize}

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

471
472
473
474
475

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

476
477

\subsubsection{Dünnbesetztheit}
478

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

481
\begin{itemize}
482
\item Ein Matrixeintrag $A_{ij} \colonequals a(\varphi_i, \varphi_j)$ ist genau dann
483
484
485
486
487
488
489
490
491
492
493
  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
494
495
496
497
498
499
500
501
502
503
504
505
506
507
\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$\\
508
     \For{alle Spalten $j$ in denen $A_{ij} \neq 0$}
509
510
511
512
513
     {
       $v_i = v_i + A_{ij} w_j$
     }
   }
  \end{algorithm}
514
Das braucht nur $O(\# \text{Nichtnulleinträge})$ Operationen (statt $O(n^2)$).
515

516
517
\bigskip

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

520
521
522
523
524
525
526
\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}
527

528
\subsection{Direkte Verfahren für dünnbesetzte Gleichungen}
529

530
531
\begin{itemize}
 \item Es gibt direkte Verfahren zum Lösen von dünnbesetzten Gleichungssystemen.
532

533
 \item Sehr kompliziert / sehr ausgefeilt.
534

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

537
538
539
540
541
542
543
544
545
546
547
548
 \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}
549
550
551
552


\subsection{Lineare iterative Verfahren}

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

556
\bigskip
557

558
559
Für \emph{lineare} iterative Verfahren~\cite{dahmen_reusken:2008} schreibt man
$Ax=b$ als Fixpunktgleichung
560
\begin{equation*}
561
 x=x+C(b-Ax)
562
\end{equation*}
563
mit einer nichtsingulären Matrix $C \in \R^{n \times n}$.
564

565
\medskip
566

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

570
571
\bigskip

572
573
Dafür machen wir jetzt eine Fixpunktiteration:
\begin{equation*}
574
  x^{k+1}\colonequals x^k+C (b-Ax^k)=(I-CA)x^k+Cb,
575
576
577
578
  \qquad
  k = 0, 1, 2, \dots
\end{equation*}

579
Der Fehler im $k$-ten Iterationsschritt ist $e^k \colonequals x^k-x^*$. Es gilt
580
581
582
\begin{equation*}
 e^{k+1}=x^{k+1}-x^*
 =
583
\underbrace{(I-CA)}_{\text{Iterationsmatrix}}e^k.
584
\end{equation*}
585
Das erklärt den Namen: Die Abbildung $e^k \mapsto e^{k+1}$ ist linear.
586
587
588
589
590

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

591
\subsubsection{Konvergenz:}
592
593

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

599
600
601
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
602
603
604
605
606
607
608
609
610
\begin{equation*}
 \abs{\lambda}\norm{v}
 =
 \norm{\lambda v}
 =
 \norm{B v}
 \le
 \norm{B}\norm{v}.
\end{equation*}
611
612

Damit erhält man:
613

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

618
619
620
621
622
623
624
625
626
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:
627
\begin{itemize}
628
 \item Das Jacobi-Verfahren
629
  \begin{equation*}
630
631
   C = D^{-1}
  \end{equation*}
632

633
634
635
636
 \item Das Gauß--Seidel-Verfahren
  \begin{equation*}
   C = (D+L)^{-1}
  \end{equation*}
637

638
 \item Das symmetrische Gauß--Seidel-Verfahren
639
  \begin{equation*}
640
   C = (D+L)^{-1}+(D+R)^{-1}-(D+R)^{-1}A(D+L)^{-1}
641
642
643
644
645
646
647
  \end{equation*}
\end{itemize}

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


\subsection{Das Jacobi-Verfahren}
648
649
650
651

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

652
Wir nehmen im Folgenden an, dass $A_{ii} \neq 0$ für alle $i=1,\ldots,n$.
653
654
655
656
657

\medskip

Betrachte das lineare Gleichungssystem:
\begin{align*}
658
659
660
    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
661
\end{align*}
662
\emph{Idee:} Löse $i$-te Zeile nach $x_i$ auf für $i=1,\ldots,n$:
663
\begin{align*}
664
665
666
    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
667
668
669
\end{align*}
Mache daraus ein iteratives Verfahren:
\begin{align*}
670
671
672
    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
673
674
675
\end{align*}
Für $i=1,\ldots,n$
\begin{equation*}
676
677
 x_i^{k+1}
 =
678
 \frac{1}{A_{ii}} \bigg(b_i-\sum_{\substack{j=1\\j\neq i}}^n A_{ij}x_j^k \bigg)
679
 =
680
 x^k_i + \frac{1}{A_{ii}} \bigg(b_i-\sum_{j=1}^n A_{ij}x_j^k \bigg).
681
\end{equation*}
682
683
684
Beachte:
\begin{itemize}
 \item Die Rechnungen für die verschiedenen $i$ sind voneinander unabhängig $\implies$ leicht zu parallelisieren.
685

686
687
 \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}
688

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

694
695
696
697
698
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}
 =
699
 x^k_i + \frac{\eta}{A_{ii}} \bigg(b_i-\sum_{j=1}^n A_{ij}x_j^k \bigg)
700
701
702
 \qquad \text{bzw.} \qquad
 x^{k+1} = x^k + \eta D^{-1} ( b - Ax^k).
\end{equation*}
703
704
705
706
707

\subsubsection{Konvergenz}

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

711
712
713
714
715
716
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}.
717

718
719
Wir greifen die Konvergenztheorie des Jacobi-Verfahrens in Kapitel~\ref{sec:jacobi_smoother}
wieder auf.
720

721
\subsubsection{Konvergenzgeschwindigkeit}
722
723
724
725
726

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

\bigskip

727
728
Sei $A$ die Matrix des Poisson-Problems auf einem uniformen Gitter
mit Dreiecks- oder Viereckselementen.  Eine Zeile von $A$ ist dann
729
730
731
732
\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*}

733
734
735
736
\begin{exercise}
 Rechnen Sie nach dass man tatsächlich diese Matrix erhält.
\end{exercise}

737

738
739
740

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

741
Wir versuchen jetzt, den Spektralradius von $I-CA = I - \eta D^{-1}A$ abzuschätzen:
742
\begin{align*}
743
744
745
746
    \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}}.
747
748
\end{align*}
Der Bruch ist gerade der \emph{Rayleigh-Quotient}
749
750
751
752
753
754
\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.

755
756
\medskip

757
Daraus folgt dass
758
759
760
761
762
763
764
765
766
767
768
769
\begin{align*}
 \rho (I - \eta D^{-1}A )
 & =
 \sup \left\lbrace \abs[\Big]{1- \frac{1}{4} \eta h^2 \lambda} \; : \; \text{$\lambda$ Eigenwert von $A$} \right\rbrace \\
 %
 & =
 1-2 \eta \sin^2 \left( \frac{1}{2} \pi h \right).
\end{align*}
Mit $\eta = 1$ erhält man:
\todo[inline]{Zeige den allgemeinen Fall!}
\begin{align*}
 \rho (I - \eta D^{-1}A )
770
 =
771
772
 \cos (\pi h).
\end{align*}
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802

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.

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

806
807
\medskip

808
Das Jacobi-Verfahren scheint also kein sonderlich gutes Verfahren zu sein.
809
810
811
812
813
814
815

\subsection{Das Gauß-Seidel-Verfahren}

Betrachte noch einmal die Jacobi-Rechenvorschrift:
\begin{equation*}
 x_i^{k+1}
 =
Sander, Oliver's avatar
Sander, Oliver committed
816
817
 \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).
818
819
820
\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
821
822
\medskip

823
824
Gauß-Seidel-Verfahren:
\begin{align*}
Sander, Oliver's avatar
Sander, Oliver committed
825
826
827
828
829
 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).
830
831
832
\end{align*}
$x_i^{k+1}$ hängt von $x_{i-1}^{k+1}$ ab $\implies$ keine Parallelisierung möglich.

833
834
835
836
837
\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}
838
839
840
841
842
843
844
845
846

\subsubsection{Konvergenz}

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

\medskip

Und in der Tat:
\begin{theorem}
847
848
849
850
Das Gauß-Seidel-Verfahren konvergiert, wenn
\begin{itemize}
 \item $A$ symmetrisch und positiv definit ist und/oder
 \item $A$ irreduzibel und diagonaldominant ist.
851
852
853
854
855
856
\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
857
 \item Es gilt $\rho(I-(D+L)^{-1}A) = ( \rho (I-D^{-1}A ) )^2$.
858
859
860
861
  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}.

862
863
864
865
866
  \begin{exercise}
   Finden Sie die genaue Stelle, und sagen Sie mir bescheid.
  \end{exercise}


867
868
869
 %
 \item Deshalb
  \begin{equation*}
Sander, Oliver's avatar
Sander, Oliver committed
870
   \rho(I-(D+L)^{-1}A) = \cos^2(\pi h)
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
  \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.

898
899
900
901
902
903
904
905
906
\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}
907
908
909
910
911
912
913
914
915

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

916
\subsubsection{Das Verfahren der konjugierten Gradienten (CG-Verfahren)}
917
918

\begin{itemize}
919
920
921
922
923
924
 \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.

925
926
927
928
 \item Geriet zwischenzeitlich in Vergessenheit
 \item Erlebte Revival als iterative Methode
\end{itemize}

929
930
\begin{theorem}[{\cite{shewchuk:1994}}]
Sei $\kappa$ die Kondition der Matrix $A$. Dann gilt nach $k$ Schritten des CG-Verfahrens
931
\begin{align*}
932
 \norm{e^k}_A \leq 2 \left( \frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1} \right)^k \cdot \norm{e^0}_A.
933
\end{align*}
934
935
\end{theorem}

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

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

941
942
 \item Aber man benötigt wieder einen Vorkonditionierer, sonst ist auch
   das CG-Verfahren zu langsam.
943

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

946
947
 \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.
948

949
950
 \item Die Konstruktion passender $C$ bildet ein großes Gebiet innerhalb
   der Numerik.
951

952
 \item Eine Möglichkeit: Mehrgitter!
953
954
\end{itemize}

955
\section{Mehrgitterverfahren: Der anschauliche Zugang}
956

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

960
\subsection{Glättung}
961

962
963
Gauß--Seidel-Verfahren und Jacobi-Verfahren sind für typische Finite-Elemente-Matrizen
extrem langsam.
964

965
966
967
968
\medskip

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

970
971
972
\smallskip

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

974
975
976
\bigskip

\emph{Beispiel:} Wir betrachten wieder das Poisson-Problem
977
978
979
980
981
\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$.
982
983
984

\medskip

985
986
987
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:
988
989
990
991
992
993
994
995
996
997
998
999
1000
\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

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

\bigskip

For faster browsing, not all history is shown. View entire blame