Numerik

Mehrgitterverfahren



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

Feingitter
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}
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 is $\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.

 

\begin{figure}\centerline{\hbox{ \psfig{file=amg_plots/freq.ps,width=10cm}} }\end{figure}
Abbildung 2:
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}$: „schlechte“ Darstellung auf dem groben Gitter und
* die Menge der niederfrequenten Gitterfunktionen $X_{h,glatt}= X_h \setminus X_{h,\rm oz}$: „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.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).

 

 

Eine Kombination aus klassischem Iterationsverfahren und der Grobgitterkorrektur $\eqref{eq6}$ ergibt das Zweigitterverfahren:

  1. Glättung: führe einige Schritte (meist weniger als vier) mit einem geeigneten Iterationsverfahren durch
  2. Grobgitterkorrektur: anschließend wird eine neue Näherung mittels $\eqref{eq5}$ berechnet
  3. falls $\| A_h x_h^{new} – f_h \| < tol$ oder eine maximale Anzahl von Iterationen überschritten ist: Stop; sonst gehe zu 1.

Algorithmus \ref{ZGMstep} beschreibt ein Unterprogramm, das einen Zweigitterschritt ausführt, und Algorithmus \ref{ZGM} die Zweigitteriteration.

\begin{algorithm}% latex2htm id marker 469\caption[Zweigitterschritt]{ zgmst......, e_H$ ; \hfill \COMMENT{ Grobgitterkorrektur }\end{algorithmic}\end{algorithm}

\begin{algorithm}% latex2htm id marker 484\caption{Zweigitterverfahren } \i......m $<$\space tol) or (step $>$\space stepmax)};\end{algorithmic}\end{algorithm}

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 $\ref{ZGMstep}$)
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 <ahref=“itnode5.htm#ZGM“>4<fontface=“Arial“> 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 ([<ahref=“itnode7.htm#Hackbusch3″>9]). Das Unterprogramm, das einen Mehrgitterschritt
ausführt kann folgendermaßen beschrieben werden

\begin{algorithm}% latex2htm id marker 516\caption[Mehrgitterschritt]{ mgmst......alpha_2$\space Nachgl\uml {a}ttungen }\ENDIF\end{algorithmic}\end{algorithm}

Die Mehrgitteriteration ist dann

\begin{algorithm}% latex2htm id marker 542\caption[Mehrgitterverfahren]{Mg-......m $<$\space tol) or (step $>$\space stepmax)};\end{algorithmic}\end{algorithm}

Die angesprochenen Themenbereiche können
in Hackbusch[9] nachgelesen werden.
Eine tiefergehende Darstellung der Mehrgitterverfahren findet
sich in Hackbusch [7].

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