Mehrgitterverfahren

Veröffentlicht am Veröffentlicht in Numerik

 



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 \mbox{ ist Knoten der Diskretisierung} \}$ ergebe ein Gleichungssystem
\begin{equation}
A_H \, x_H = f_H . \label{eq2}
\end{equation}

mit positiv definiter Matrix $A_H$. Durch Halbierung der Schrittweite $H$, d.h. Verfeinerung des Gitters $\Omega_H$, erhält man ein weiteres Gleichungssystem, das die gleiche Struktur wie $\eqref{eq1}$, jedoch die 4-fache Dimension, aufweist. Die Approximation der Lösung von $\eqref{eq1}$ durch den Ansatz verbessert sich, jedoch wächst der Aufwand zur Lösung des entsprechenden Gleichungssystems.

Grobgitter
Verfeinerung $h=\frac{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}
A_h \, e_h = d_h, \qquad \mbox{ bzw. } \qquad e_h = A_h^{-1} d_h . \label{eq3}
\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}
Im Vergleich zur Ausgangssituation ist nichts gewonnen, da die Lösung von $\eqref{eq3}$ genauso schwierig (aufwendig) ist, wie die der Gleichung $A_h \, x_h = f_h$. 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}
und nach $e_H$ aufzulösen: $e_H = A_H^{-1}r \, d_h$. 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}
wobei die Abbildung $p$ (die Prolongation) die Korrektur $e_H$ auf das feine Gitter überträgt.
Sind $r_h$ und $p_h$ die Matrixdarstellungen der Restriktion und Prolongation, so wird häufig
\begin{equation}
r_h = \sigma \, p_h^T \label{eq7}
\end{equation}
gewählt, dabei ist $\sigma$ ein Skalierungsfaktor. 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.

Eindimensionales Beispiel: (oben) glatte (niederfrequente) Gitterfunktion, gut darstellbar auf dem groben Gitter; (unten) hochfrequente Gitterfunktion, nicht “sichtbar” im groben Gitter

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 Menge der hochfrequenten Gitterfunktionen

$$ X_{h,{\rm oz}} \; : \quad \text{“schlechte” Darstellung auf dem groben Gitter} $$

und

die Menge der niederfrequenten Gitterfunktionen

$$X_{h,glatt}= X_h \setminus X_{h,{\rm oz}} \; : \quad \text{“gute” Darstellung auf dem groben Gitter} $$

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}

Da es (hochfrequente) Feingitterfunktionen gibt, die das grobe Gitter nicht “sieht”, bildet das durch $\eqref{eq6}$ gegebene Verfahren i. Allg. 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ämpft, die niederfrequenten Anteile in der Lösung hingegen nur langsam.
Fehler $e_h(x)$, $y$ ist die Nummer der Jacobi-Iteration.

 

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: Stopp; sonst gehe zu 1.

Algorithmus 3 beschreibt ein Unterprogramm, das einen Zweigitterschritt ausführt, und Algorithmus 4 die Zweigitteriteration.

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}

Ausgehend von einer Näherung $x_i$ kann ein Schritt (also die Berechnung der nächsten Iterierten) im eben beschriebenen Zweigitterverfahren (Algorithmus 4) 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 $\cal{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}

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 ist.

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 tiefer gehende 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

1.
Hackbusch W (1985) Multi-Grid Methods and Applications. Springer Series in Computational Mathematics. doi: 10.1007/978-3-662-02427-0
2.
Hackbusch W (1993) Iterative Lösung großer schwachbesetzter Gleichungssysteme. doi: 10.1007/978-3-663-05633-1