Mehrgitterverfahren sind in der Lage große Gleichungssysteme mit mehreren Millionen Unbekannten, die sich aus der Diskretisierung von physikalischen Simulationsmodellen ergeben und bestimmte Anforderungen erfüllen, schnell zu lösen. Im Gegensatz zu klassischen Iterationsverfahren, deren Konvergenzrate mit zunehmender Feinheit des Gitters abnimmt, behalten Mehrgittermethoden ihre guten Konvergenzeigenschaften auch für sehr feine Gitter bei.
Die Diskretisierung der Randwertaufgabe
\begin{equation}
– \Delta u ({\bf x}) = f ({\bf x}) , \qquad {\bf x}\in \Omega , \qquad u\Big|_{\partial \Omega} = g \label{eq1}
\end{equation}
auf einem Gitter \(\Omega_{H} :=\{ {\bf x}^i \; :\; {\bf x}^i \text{ ist Knoten der geometrischen Modellierung} \}\) führt beispielsweise mittels FEM auf ein Gleichungssystem der Form:
\begin{equation} \label{eq2}
A_H \, x_H = f_H .
\end{equation}
Die Matrix \(A_H\) ist positiv definit. Wird die Gitterweite \(H\) halbiert, so ergibt sich ein verfeinertes Gitter \(\Omega_h\). Das daraus resultierende Gleichungssystem besitzt dieselbe Struktur wie Gleichung \(\eqref{eq2}\), jedoch mit deutlich höherer Dimension (typischerweise etwa viermal so groß). Zwar verbessert sich die Näherung der Lösung durch diese Verfeinerung, jedoch steigt auch der Rechenaufwand deutlich an.
Verfeinerung mit $h=H/2$
Sei nun $x_h^{\star}$ die Lösung des Gleichungssystems $A_h \, x_h = f_h$ im Gitter $\Omega_h$ und $\tilde{x}_h$ eine Näherung der
Lösung. Es gilt somit
$$ A_h (\tilde{x}_h – x_h^{\star}) = A_h \, \tilde{x}_h – f_h := d_h.$$
Der Vektor $d_h$ ist der Defekt und
$$
e_h := \tilde{x}_h – x_h^{\star}
$$
der Fehler zwischen Näherung und Lösung. Insbesondere folgt die Defektgleichung
\begin{equation} \label{eq3}
A_h \, e_h = d_h, \qquad \mbox{ bzw. } \qquad e_h = A_h^{-1} d_h .
\end{equation}
Ist durch Lösung des Ersatzgleichungssystems der Fehler \(e_h\) bekannt, so ergibt sich die gesuchte Lösung sofort durch
\begin{equation}
x_h^{\star} = \tilde{x}_h – e_h. \label{eq4}
\end{equation}
Restriktion und Prolongation
Die (Defekt-) Gleichung \(\eqref{eq3}\) ist im Prinzip genauso schwer zu lösen wie das ursprüngliche Gleichungssystem \(A_h \, x_h = f_h\). Es wäre also nichts gewonnen. Der entscheidende Gedanke ist nun, die Defektgleichung auf das grobe Gitter $\Omega_H$ mittels einer Abbildung $r: \Omega_h \longrightarrow \Omega_H$ (Restriktion) zu übertragen:
\begin{equation}
A_H \, e_H = (r \, d_h) \label{eq5} .
\end{equation}
Die (Grob-) Lösung lautet dann : $e_H = A_H^{-1} \, r \, d_h$. Dieser korrigierende Anteil \(e_H\) wird anschließend auf das feine Gitter zurückprojiziert.
Der Defekt $d_h$ wird also (in geeigneter Weise) auf das grobe Gitter „transportiert“. Die anschließende Lösung des Gleichungsystems $\eqref{eq5}$ ist weniger aufwendig, als die der ursprünglichen Gleichung. Für ein 2D-Problem hat man $\dim A_H = (\dim A_h)/4.$
Eine neue Näherung ergibt sich nun durch
\begin{equation}
x_h^{\rm new} := \tilde{x}_h – p \, e_H . \label{eq6}
\end{equation}
Dabei ist \(p\) die sogenannte Prolongation, also die Abbildung, die die Korrektur \(e_H\) vom Grobgitter zurück auf das feine Gitter überträgt.
Sind \(r_h\) und \(p_h\) die Matrixdarstellungen von Restriktion und Prolongation, so wählt man oft:
\[
r_h = \sigma \, p_h^T \tag{7} \label{eq7} .
\]
Hierbei ist \(\sigma\) ein geeigneter Skalierungsfaktor, der die Konsistenz und Stabilität der Methode verbessert. Implizit geht hier die Wahl der Diskretisierung (Finite-Elemente-Methode, Finite-Differenzen-Verfahren) und der Skalarprodukte
für die Vektoren auf den verschiedenen Gittern ein.
 |
Gitterfunktionen
Ausgehend vom groben Gitter $\Omega_H$ läßt sich die Menge der Vektoren $X_h$ (auch Menge der Gitterfunktionen genannt) zum feinen Gitter $\Omega_h$ in zwei Kategorien einordnen
- die Hochfrequente Gitterfunktionen \(X_{h,\mathrm{oz}}\): $ \;$ Diese Anteile lassen sich auf dem groben Gitter schlecht darstellen.
- die Niederfrequente Gitterfunktionen \(X_{h,\mathrm{glatt}} := X_h \setminus X_{h,\mathrm{oz}}\): $ \;$ Diese Anteile sind auf dem groben Gitter gut aufgelöst.
Ein Vektor $x_h \in X_h$ kann also in hoch- und niederfrequente Anteile zerlegt werden
\begin{equation}
x_h = \sum_{e_i \in X_{h,\rm glatt}} \alpha_i \, e_i + \sum_{e_j \in X_{h.\rm oz}} \alpha_j \, e_j .
\end{equation}
Glättung
Da es (hochfrequente) Feingitterfunktionen gibt, die das grobe Gitter nicht „sieht“, bildet das durch $\eqref{eq6}$ gegebene Verfahren i.a. keine konvergente Iteration. Die folgende Beobachtung hilft an dieser Stelle entscheidend weiter: Klassische Iterationsverfahren, wie etwa das Gauß-Seidel-Verfahren, haben gerade Probleme die niederfrequenten Anteile in der Lösung zu erfassen, reduzieren aber die hochfrequenten Fehleranteile gut.
Darstellung des Fehler zu einem 1D-Modellproblem bgzl. des Jacobi-Verfahrens : der hochfrequenten Anteil wird schnell gedämft, die niederfrequenten Anteile in der Lösung hingegen nur langsam (unten: $e_h(x,y)$, wobei $y$ die Nummer der Jacobi-Iteration ist).
Zweigitterverfahren
Eine Kombination aus klassischem Iterationsverfahren und der Grobgitterkorrektur $\eqref{eq6}$ ergibt das Zweigitterverfahren:
- Glättung: führe einige Schritte (meist weniger als vier) mit einem geeigneten Iterationsverfahren durch
- Grobgitterkorrektur: anschließend wird eine neue Näherung mittels $\eqref{eq5}$ berechnet
- falls $\| A_h x_h^{new} – f_h \| < tol$ oder eine maximale Anzahl von Iterationen überschritten ist: Stop; sonst gehe zu 1.
Algorithmus 3 beschreibt ein Unterprogramm, das einen Zweigitterschritt ausführt
\[
\begin{array}{ll}
\textbf{Algorithmus 3: zgmstep}(x_h, f_h, l) & \text{// Zweigitterschritt mit l Glättungen} \\
x_h = S_l(x_h, f_h) & \text{// Glättung (l Schritte)} \\
d_h = A_h \cdot x_h – f_h & \text{// Defekt} \\
d_H = r \cdot d_h & \text{// Restriktion auf das Grobgitter} \\
e_H = A_H^{-1} \cdot d_H & \text{// Lösung auf dem Grobgitter} \\
x_h = x_h – p_h \cdot e_H & \text{// Grobgitterkorrektur} \\
\end{array}
\]
und Algorithmus $4$ die Zweigitteriteration
\[
\begin{array}{ll}
\textbf{Algorithmus 4: Zweigitterverfahren} & \\
\text{Initialisierung …} & \\
\textbf{repeat} & \\
\quad \text{zgmstep}(x_h, f_h, l) & \text{// ein Zweigitterschritt} \\
\quad \text{norm} = \| A_h \cdot x_h – f_h \| & \text{// Fehlerberechnung} \\
\quad \text{step} = \text{step} + 1 & \\
\textbf{until} (\text{norm} < \text{tol}) \; \text{oder} \; (\text{step} > \text{maxsteps}) &
\end{array}
\]
Dabei ist
\begin{equation}
{\cal S}_l (x_h,f_h) := S_h \ x_h + N_h \, f_h
\end{equation}
eine konsistente Iteration, die zur Glättung benutzt wird. Ist diese die Jacobi-Iteration, so gilt
\begin{equation}
S_h =I_h – D_h^{-1} A_h , \qquad N_h=D_h^{-1}, \quad D_h:=diag(A_h).
\end{equation}
Iterationsmatrix
Ausgehend von einer Näherung $x_i$ kann ein Schritt (also die Berechnung der nächsten Iterierten) im eben beschriebenen Zweigitterverfahren auch als lineares Iterationsverfahren
\begin{equation}
x^{(i+1)} = S_h x^{(i)} – p_h A_H^{-1} r_h \left(A_h S_h x_h^{(i)} – \tilde{f}_h \right)
\end{equation}
geschrieben werden. Die Iterationsmatrix des Zweigitterverfahrens ist also
\begin{equation}
M_h^{\rm ZGM} := \left( I_h – p_h \ A_H^{-1} \, r_h \, A_h \right) \, S_h.
\end{equation}
Mit Hilfe von Annahmen über die Regularität der Lösung der zu Grunde liegenden partiellen Differentialgleichung, weiterer Voraussetzungen an die Diskretisierung (geschachtelte Unterräume, Approximationseigenschaften) und an die Glättungsiteration $\mathcal{S}_l$ kann gezeigt werden, daß der Spektralradius der Iterationsmatrix $M_h^{\rm ZGM}$ unabhängig von der Schrittweite h ist:
\begin{equation}
\rho( M_h^{\rm ZGM}) \leq \zeta(\alpha) < 1 , \qquad \mbox{f.a.} \quad \alpha \geq \alpha_0.
\end{equation}
Mehrgitterverfahren
Das Verfahren 4 konvergiert demnach, unter hier nicht näher genannten Voraussetzungen, unabhängig von $h$. Der Übergang vom Zweigitter- zum Mehrgitterverfahren ist jetzt nur noch ein kleiner Schritt: das exakte Lösen der Grobgittergleichung wird nun wiederum durch ein Zweigitterverfahren ersetzt, usw. Die sehr guten Konvergenzeigenschaften bleiben dabei erhalten und es zeigt sich, daß der Aufwand optimal. Das Unterprogramm, das einen Mehrgitterschritt ausführt kann folgendermaßen beschrieben werden
Die Mehrgitteriteration ist dann
Die angesprochenen Themenbereiche können in Hackbusch[2] nachgelesen werden. Eine tiefergehende Darstellung der Mehrgitterverfahren findet sich in Hackbusch [1].
Beispiel: 2D-Modellproblem mit $\Omega=[0,1]$, $g=0$ und $f(x,y)= sin(x \, \pi)* sin(y*\pi)$
- Diskretisierung: 33×33 Gitter
- Startvektor: x=random
- maxlevel=5, $\gamma=1$
- Glättung: zwei Gauß-Seidel-Iterationen
- Konvergenzrate: 0.2
Literatur
Iterative Lösung großer schwachbesetzter Gleichungssysteme Buch
Teubner, Stuttgart, 1991.
Multi-grid methods and applications Buch
Springer, Berlin, 1985, ISBN: 978-3-540-12761-1.