4.6. Gesamt- und Einzelschrittverfahren

Viele mathematische Probleme, die numerisch gelöst werden sollen, führen zu einem großen linearen Gleichungssystem \(Au = f\), bei dem die Matrix nur an relativ wenigen Stellen Einträge ungleich Null besitzt. Dies tritt vor allem bei der Diskretisierung von Differentialgleichungen (DGL) auf, wie das folgende motivierende Beispiel demonstriert.

Beispiel 4.5 (Diskretisierung einer Helmholtz-artigen DGL)

Wir betrachten ein eindimensionales Randwertproblem auf einem Intervall \(\Omega = [0,1]\) für \(\sigma > 0\) der Form:

(4.8)\[\begin{split} - u''(x) + \sigma u(x) \ &= \ f(x), \qquad x \in (0,1), \\ u(0) \ = \ u(1) \ &= \ 0. \qquad \text{(Dirichlet-Randbedingungen)} \end{split}\]

Diese Helmholtz-artige Differentialgleichung wird beispielsweise zur Entrauschung von Daten in der mathematischen Bildverarbeitung genutzt.

Die kontinuierliche Differentialgleichung in (4.8) lässt sich im Eindimensionalen auf einem diskreten Gitter

\[\Omega_h \ \coloneqq \ \lbrace x_i \in [0,1] \: | \: x_i = i \cdot h \rbrace\]

mit Schrittweite \(h \coloneqq 1 / N\) und \(N+1 \in \N\) Gitterpunkten approximieren (siehe Vorlesung Numerik 2: Diskretisierung und Optimierung). Dies führt zu einem linearen Gleichungssystem mit der (bis auf den Rand) unbekannten Lösung \(u \in \R^{N+1}\) und der Diskretisierung \(f_i \coloneqq f(x_i)\) der Form

(4.9)\[\begin{split} \frac{-u_{i-1} + 2u_i - u_{i+1}}{h^2} + \sigma u_i \ &= \ f_i, \qquad i=1,\ldots,N-1, \\ u_0 \ = \ u_N \ &= \ 0. \end{split}\]

Dies können wir kompakt als lineares Gleichungssystem der Form \(Au = f\), wobei nun \(A \in \R^{(N-1) \times (N-1)}\) eine Tridiagonalmatrix und \(u,f \in \R^{N-1}\) ist mit

\[A \ \coloneqq \ \frac{1}{h^2} \begin{pmatrix} 2 + \sigma h^2 & -1 & 0 & \cdots & 0\\ -1 & 2 + \sigma h^2 & -1 & \ddots & \vdots\\ 0 & -1 & \ddots & \ddots & \vdots&\\ \vdots & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & 0 & -1 & 2 + \sigma h^2 \end{pmatrix}, \ u \ \coloneqq \ \begin{pmatrix} u_1\\ \vdots\\ u_{N-1} \end{pmatrix}, \ f \ \coloneqq \ \begin{pmatrix} f_1\\ \vdots\\ f_{N-1} \end{pmatrix}.\]

Eine Lösung \(u \in \R^{N-1}\) des linearen Gleichungssystems (4.9) stellt dann eine diskrete Approximation einer Lösung der Differentialgleichung (4.8) dar.

Die Matrix \(A\) in Beispiel 4.5 ist von einer besonderen Gestalt, da sie symmetrisch ist und zum größten Teil nur Nullen als Einträge hat. Solche Matrizen nennt man auch dünn-besetzt (im Englischen: sparse). Man kann außerdem zeigen, dass die Matrix positiv-definit ist für \(\sigma > 0\). Nun könnte man auf Grund dieser Eigenschaften eine direkte Lösung mittels der Cholesky-Zerlegung aus Cholesky-Verfahren umsetzen. Dies ist in der Regel jedoch weniger effizient für dünn-besetzte Matrizen, für die sich häufig iterative Lösungsverfahren besser eignen. Solche Verfahren wollen wir uns im Folgenden genauer anschauen.

4.6.1. Gesamtschrittverfahren

Das Gesamtschrittverfahren, das oft auch Jacobi-Verfahren genannt wird, basiert auf einer relativ einfachen Idee. Möchte man nämlich ein allgemeines lineares Gleichungssystem

\[Au \ = \ f, \qquad A \in \R^{n \times n}, \quad u,f \in \R^n,\]

lösen, wobei \(A\) eine dünn-besetzte Matrix darstellt, so zerlegt man die Matrix (additiv) in folgende Matrizen:

\[A \ = \ L + D + R, \qquad L,D,R \in \R^{n \times n}.\]

Hierbei stellt \(D\) die Hauptdiagonale der Matrix \(A\) dar, und \(L\) und \(R\) sind jeweils linke untere und rechte obere Dreiecksmatrizen, die sich aus den entsprechenden Einträgen aus \(A\) ergeben.

Damit lässt sich das allgemeine lineare Gleichungssystem nun wie folgt umschreiben:

\[Au \ = \ (L + D + R)u \ = \ Du + (L+R)u \ = \ f.\]

Wir approximieren die obige Gleichung nun durch ein iteratives Verfahren der Form

\[D u^{k+1} + (L + R) u^k \ \approx \ f,\]

für einen geeigneten Startpunkt \(x_0 \in \R^n\). Falls für die Einträge der Diagonalmatrix \(D\) gilt \(d_{i,i} \neq 0, i=1,\ldots,n\) ist eine Inversion der Matrix \(D\) wohldefiniert und darüber hinaus sehr einfach zu realisieren und wir erhalten

(4.10)\[u^{k+1} \ = \ D^{-1}(f - (L + R) u^k ).\]

Wenn man diese Iterationsvorschrift für die einzelnen Einträge des Lösungsvektors \(u^{k+1} \in \R^n\) betrachtet ergibt sich folgende kompakte Form:

(4.11)\[u_i^{k+1} \ = \ \frac{1}{a_{i,i}}\left(f_i - \sum_{j\neq i} a_{i,j} u_j^k \right), \qquad i=1,\ldots,n.\]

Die Berechnung jedes Eintrags ist offensichtlich unabhängig von den anderen Einträgen und kann somit effizient parallelisiert werden. Im Gesamtschrittverfahren werden also pro Iteration erst alle Einträge des Lösungsvektors \(u\) aktualisiert und anschließend für die nächste Iteration verwendet.

Bemerkung 4.8 (Rechenaufwand des Gesamtschrittverfahrens)

Betrachtet man die Rechenvorschrift des Gesamtschrittverfahrens in (4.11), so ist klar, dass man im schlechtesten Fall einer voll-besetzten Matrix \(A\) für die Berechnung eines Eintrags \(u_i^{k+1}\) genau \(n\) FLOPS benötigt und somit für die Berechnung der Iterierten \(u^{k+1} \in \R^n\) insgesamt \(n^2\) FLOPS benötigt.

Nehmen wir an, dass wir die echte Lösung \(u \in \R^n\) des linearen Gleichungssystems in \(m \in \N\) Iterationsschritten approximieren wollen, so ist der Rechenaufwand für das Gesamtschrittverfahren \(m \cdot n^2\) und liegt somit in \(\mathcal{O}(n^2)\). Für wenige Iterationen \(m << n\) ist dies schon deutlich effizienter als direkte Lösungsverfahren, die wir in Direkte Lösung linearer Gleichungssysteme diskutiert haben.

Dazu kommt noch der Vorteil, dass für dünn-besetzte Matrizen mit nur \(k \in \N\), \(k << n\) Einträgen ungleich Null pro Zeile der Rechenaufwand nochmal signifikant sinkt und man insgesamt \(m \cdot k \cdot n\) Rechenoperationen zur Approximation einer Lösung \(\overline{u} \in \R^n\) benötigt.

Beispiel 4.6 (Gesamtschrittverfahren für Helmholtz-artige DGL)

Wenden wir das Gesamtschrittverfahren zur Lösung des linearen Gleichungssystems für die Helmholtz-artige Differentialgleichung aus Beispiel 4.5 an, so erhalten wir die folgende Iterationsvorschrift für die unbekannte Lösung \(u \in \R^{N-1}\):

\[\begin{split} u_1^{k+1} \ &= \ \frac{1}{2 / h^2 + \sigma}\left(f_i +\frac{1}{h^2}u_{2}^k \right), \\ u_i^{k+1} \ &= \ \frac{1}{2 / h^2 + \sigma}\left(f_i +\frac{1}{h^2}(u_{i-1}^k + u_{i+1}^k) \right), \qquad i=2,\ldots,N-2,\\ u_{N-1}^{k+1} \ &= \ \frac{1}{2 / h^2 + \sigma}\left(f_i +\frac{1}{h^2}u_{N-2}^k \right).\ \end{split}\]

Da wir nur maximal \(3\) Einträge ungleich Null in der Matrix \(A\) haben, lässt sich eine Approximation der Lösung \(u\) in weniger als \(3 \cdot m \cdot N\) FLOPS berechnen, wobei \(m \in \N\) die Anzahl der Iterationen darstellt.

Nun stellt sich die Frage wann das Gesamtschrittverfahren gegen eine Lösung \(\overline{u} \in \R^n\) des linearen Gleichungssystems \(Au = f\) konvergiert. Wir werden zeigen, dass dies durch eine spezielle Klasse von Matrizen, die wir im Folgenden einführen werden, immer der Fall ist.

Definition 4.4 (Diagonaldominante Matrizen)

Eine Matrix \(A \in \R^{n \times n}\) heißt stark diagonaldominant, falls die Beträge ihrer Diagonalelemente \(a_{i,i}, i=1,\ldots,n\) jeweils größer sind als die Summer der Beträge der anderen jeweiligen Zeileneinträge \(a_{i,j}\), d.h. wenn gilt

\[\sum_{j=1, j\neq i}^n |a_{i,j}| \ < \ |a_{i,i}|, \qquad \forall \: i=1,\ldots,n.\]

Man nennt diese Bedingung auch das starke Zeilensummenkriterium.

Ist die Ungleichung nicht strikt, d.h. es gilt lediglich

\[\sum_{j=1, j\neq i}^n |a_{i,j}| \ \leq \ |a_{i,i}|, \qquad \forall \: i=1,\ldots,n.\]

so ist die Matrix schwach diagonaldominant.

Bemerkung 4.9 (Spaltensummenkriterium)

Analog zum Zeilensummenkriterium in Definition 4.4 gibt es auch ein starkes und schwaches Spaltensummenkriterium. Dieses ist offensichtlich äquivalent zum Zeilensummenkriterium der transponierten Matrix.

Für stark diagonaldominante Matrizen können wir nun die Konvergenz des Gesamtschrittverfahrens nachweisen.

Satz 4.5 (Konvergenz des Gesamtschrittverfahrens)

Das Gesamtschrittverfahren zur Lösung des linearen Gleichungssystems \(Au=f\) für eine reguläre Matrix \(A\in \R^{n\times n}\) und \(u,f \in \R^n\) konvergiert für jeden Startvektor \(u^0 \in \R^n\) gegen eine Lösung, falls \(A\) das starke Zeilensummenkriterium oder das starke Spaltensummenkriterium erfüllt.

Proof. Wir verwenden die additive Zerlegung der Matrix \(A\) in \(A = L + D +R\) und sehen ein, dass \(D\) nur Einträge ungleich Null auf der Hauptdiagonale besitzt, da \(A\) sonst nicht diagonaldominant wäre.

Um die Konvergenz des Gesamtschrittverfahrens zu zeigen, müssen wir zeigen, dass die Fixpunktiteration

\[u^{k+1} \ = \ G(x^k) \ \coloneqq \ -D^{-1}(L+R)u^k + D^{-1}f\]

kontraktiv ist. Hierzu definieren wir uns die Matrix \(B \in \R^{n \times n}\) als \(B \coloneqq - D^{-1}(L+R)\). Aus der Gestalt von \(B\) wird klar, dass für die Einträge \(B_{i,j} \in \R\) gilt:

\[b_{i,j} \ = \ \begin{cases} \ 0, \quad &\text{ falls } i = j,\\ \ \frac{a_{i,j}}{a_{i,i}}, \quad &\text{ sonst. } \end{cases}\]

Falls \(A\) das starke Zeilensummenkriterium erfüllt, folgt für die Matrix \(B\) schon direkt

\[\begin{split} ||B||_\infty \ &= \ \max_{i=1,\ldots,n} | \sum_{j=1}^n b_{i,j} | \ = \ \max_{i=1,\ldots,n} | \sum_{j \neq i} \frac{a_{i,j}}{a_{i,i}} | \\ &\leq \ \max_{i=1,\ldots,n} \frac{ \sum_{j \neq i} |a_{i,j}|}{|a_{i,i}|} \ < \ \max_{i=1,\ldots,n} \frac{ |a_{i,i}|}{|a_{i,i}|} \ = \ 1. \end{split}\]

Analog kann man zeigen, dass für die Spaltensummernnorm \(|| B ||_1 < 1\) gilt, falls \(A\) das starke Spaltensummenkriterium erfüllt.

Nun können wir die Kontraktivität der Fixpunktiteration \(G\) leicht zeigen, denn es gilt:

\[\begin{split} || G(u^k) - G(\overline{u})|| \ &= \ || \underbrace{- D^{-1}(L+R)}_{= \, B} u^k + D^{-1}b - (\underbrace{- D^{-1}(L+R)}_{= \, B}\overline{u} + D^{-1}b)|| \\ &= \ || Bu^k - B\overline{u}|| \ \leq \ \!\! \underbrace{||B||}_{=: \, q \, < \, 1} \!\!\! \cdot ||u^k - \overline{u}||. \end{split}\]

Nun können wir den Banachschen Fixpunktsatz aus Satz 4.1 anwenden, der uns die Konvergenz des Gesamtschrittverfahrens liefert. ◻

4.6.2. Einzelschrittverfahren

Wenn man auf eine parallele Berechnung der Einträge \(u_i^{k+1}, i=1,\ldots,n\) in (4.11) verzichtet und statt dessen eine sequentielle Berechnung des Vektors \(u^{k+1} \in \R^n\) vornimmt, so ergibt sich hieraus eine weitere Möglichkeit. Anstatt die Einträge vollkommen unabhängig voneinander zu berechnen könnte man auch auf die bereits aktualisierten, vorangehenden Einträge zurückgreifen, da diese schließlich der echten Lösung \(u \in \mathbb{R}^n\) näher sein sollten.

Dies ist die grundlegende Idee beim Einzelschrittverfahren, das auch unter dem Namen Gauss–Seidel-Verfahren bekannt ist. Analog zum Gesamtschrittverfahren in Gesamtschrittverfahren zerlegen wir die gegebene Matrix \(A\) wieder additiv zu \(A = L + D + R\), jedoch schreiben wir das lineare Gleichungssystem nun wie folgt auf:

\[Au \ = \ (L + D + R)u \ = \ (L + D)u + Ru \ = \ f.\]

Um nun also auf die bereits aktualisierten Einträge der approximativen Lösung \(u^{k+1} \in \R^n\) zurück zu greifen, approximieren wir das obige lineare Gleichungssystem durch

(4.12)\[(L + D)u^{k+1} + Ru^k \ \approx \ f.\]

Die Matrix \((L+D) \in \R^{n\times n}\) ist nun eine linke untere Dreiecksmatrix, die genau dann invertierbar ist, falls für die Einträge der Diagonalmatrix \(D\) gilt \(d_{i,i} \neq 0, i=1,\ldots,n\). Unter dieser Annahme lässt sich dann eine Iterationsvorschrift für das Einzelschrittverfahren in Matrixschreibweise angeben als

\[x^{k+1} \ = \ (L + D)^{-1}(f - Ru^k).\]

Wenn man diese Iterationsvorschrift für die einzelnen Einträge des Lösungsvektors \(u^{k+1} \in \R^n\) betrachtet ergibt sich aus (4.12) für alle \(i=1,\ldots,n\) folgende Form:

(4.13)\[\begin{split} & \sum_{j=1}^{i-1} a_{i,j} u_j^{k+1} + a_{i,i} u_i^{k+1} + \sum_{j=i+1}^{n} a_{i,j} u_j^{k} \ = \ f \\ \Rightarrow \quad & a_{i,i} u_i^{k+1} \ = \ f - \sum_{j=1}^{i-1} a_{i,j} u_j^{k+1} - \sum_{j=i+1}^{n} a_{i,j} u_j^{k}\\ \Rightarrow \quad & u_i^{k+1} \ = \ \frac{1}{a_{i,i}}\left(f_i - \sum_{j=0}^{i-1} a_{i,j} u_j^{k+1} - \sum_{j = i+1}^n a_{i,j} u_j^k \right). \end{split}\]

Bemerkung 4.10 (Rechenaufwand des Einzelschrittverfahrens)

Betrachtet man die Rechenvorschrift des Gesamtschrittverfahrens in (4.13), so ist klar, dass man analog zum Gesamtschrittverfahren im schlechtesten Fall einer voll-besetzten Matrix \(A\) für die Berechnung eines Eintrags \(u_i^{k+1}\) genau \(n\) FLOPS benötigt und somit für die Berechnung der Iterierten \(u^{k+1} \in \R^n\) insgesamt \(n^2\) FLOPS benötigt. Der Rechenaufwand beträgt also genauso \(m\cdot n^2\) FLOPS für \(m\in N\) Iterationsschritte wie beim Gesamtschrittverfahren in Bemerkung 4.8, lässt sich allerdings nicht mehr parallelisieren.

Beispiel 4.7 (Einzelschrittverfahren für Helmholtz-artige DGL)

Wenden wir das Einzelschrittverfahren zur Lösung des linearen Gleichungssystems für die Helmholtz-artige Differentialgleichung aus Beispiel 4.5 an, so erhalten wir die folgende Iterationsvorschrift für die unbekannte Lösung \(u \in \R^{N-1}\):

\[\begin{split} u_1^{k+1} \ &= \ \frac{1}{2 / h^2 + \sigma}\left(f_i +\frac{1}{h^2}u_{2}^k \right), \\ u_i^{k+1} \ &= \ \frac{1}{2 / h^2 + \sigma}\left(f_i +\frac{1}{h^2}(u_{i-1}^{k+1} + u_{i+1}^k) \right), \qquad i=2,\ldots,N-2,\\ u_{N-1}^{k+1} \ &= \ \frac{1}{2 / h^2 + \sigma}\left(f_i +\frac{1}{h^2}u_{N-2}^{k+1} \right).\ \end{split}\]

Da wir nur maximal \(3\) Einträge ungleich Null in der Matrix \(A\) haben, lässt sich eine Approximation der Lösung \(u\) in weniger als \(3 \cdot m \cdot N\) FLOPS berechnen, wobei \(m \in \N\) die Anzahl der Iterationen darstellt.

Bemerkung 4.11 (Konvergenz von Gesamt- und Einzelschrittverfahren)

Sowohl das Gesamt- als auch das Einzelschrittverfahren Verfahren konvergieren linear gemäß Definition 4.2. Dennoch bevorzugt man häufig das Einzelschrittverfahren, da es insgesamt schneller konvergiert als das Gesamtschrittverfahren und man somit Iterationsschritte einspart. Die Konvergenzgeschwindigkeit des Einzelschrittverfahrens lässt sich sogar noch weiter verbessern indem man ein sogenanntes Successive-Over-Relaxation (SOR) Verfahren verwendet.

Die Diagonaldominanz der Matrix \(A\) des linearen Gleichungssystems \(Au=f\) ist ebenfalls eine hinreichendes Kriterium für die Konvergenz des Einzelschrittverfahrens, analog zur Aussage des Gesamtschrittverfahrens in Satz 4.5.