Gauß-Seidel-Verfahren

Grundidee und Ablauf

Das Gauß-Seidel-Verfahren ist ein iteratives Verfahren zur näherungsweisen Lösung eines linearen Gleichungssystems mit n Variablen. Anders als ein exaktes Lösungsverfahren erzeugt es aus einem Startvektor schrittweise bessere Näherungen für den Lösungsvektor.

Für die Gleichungen a₁₁x₁+…+a₁ₙxₙ=b₁, …, aₙ₁x₁+…+aₙₙxₙ=bₙ wird im Iterationsschritt m+1 jede Gleichung nach ihrer zugehörigen Variablen aufgelöst:

x_k^(m+1) = 1/a_kk · (b_k − Σ_(i=1)^(k−1) a_ki x_i^(m+1) − Σ_(i=k+1)^n a_ki x_i^(m)), k=1,…,n.

Entscheidend ist: Bereits berechnete Werte des aktuellen Durchgangs x_i^(m+1) werden sofort verwendet. Für noch nicht berechnete Variablen werden die Werte des vorherigen Durchgangs x_i^(m) benutzt. Beim Jacobi-Verfahren würden dagegen ausschließlich Werte des vorherigen Iterationsschritts eingesetzt.

So entsteht die Folge x^(0), x^(1), x^(2), … → x. Sie nähert sich dem Lösungsvektor x nur dann an, wenn das Verfahren konvergiert. Alle Diagonalelemente a_kk müssen von 0 verschieden sein, weil durch sie dividiert wird.

Abbruch und Rechenaufwand

Der Startvektor x^(0) kann willkürlich gewählt werden. In jedem Durchgang wird zusätzlich der größte Betrag einer Änderung gespeichert:

fehler := max(fehler, |x_k^(m+1) − x_k^(m)|).

Nach dem Berechnen aller n Komponenten wird m um 1 erhöht. Das Verfahren endet, wenn fehler < fehlerschranke gilt. Die Fehlerschranke beschreibt also, wie klein die letzte Änderung des Variablenvektors sein muss.

Das Verfahren ist inhärent sequentiell: Bevor eine Gleichung bearbeitet werden kann, müssen die Ergebnisse der vorherigen Gleichungen vorliegen. Daher eignet es sich nur eingeschränkt für Parallelrechner, insbesondere bei dünnbesetzten Matrizen. Bei solchen Matrizen, die viele Nulleinträge besitzen, sinkt jedoch der Aufwand pro Iteration deutlich.

Matrixform

Für das Gleichungssystem A·x=b wird A in drei Teile zerlegt:

A = L + D + U.

Dabei ist D die Diagonalmatrix, L eine strikte untere Dreiecksmatrix und U eine strikte obere Dreiecksmatrix. Ein Iterationsschritt erfüllt zunächst

D x^(neu) = b − L x^(neu) − U x^(alt).

Daraus folgt

(D+L)x^(neu) = b − Ux^(alt)

und formal

x^(neu) = (D+L)⁻¹(b−Ux^(alt)).

Mit einem Startvektor x^(0) lautet die Iterationsvorschrift:

x^(k+1) = −(D+L)⁻¹U x^(k) + (D+L)⁻¹b.

Der Ausdruck −(D+L)⁻¹U heißt Iterationsmatrix.

Beispiel

Gegeben sind

A = ((16, 3), (7, −11)) und b = (11, 13)ᵀ.

Die Zerlegung in eine untere Dreiecksmatrix einschließlich Diagonale und eine strikte obere Dreiecksmatrix ist

A = L_* + U = ((16, 0), (7, −11)) + ((0, 3), (0, 0)).

Mit T = −L_⁻¹U und C = L_⁻¹b wird x^(k+1)=Tx^(k)+C. Näherungsweise ergeben sich

T = ((0,000, −0,1875), (0,000, −0,1194))

und C = (0,6875, −0,7439)ᵀ.

Für den Startvektor x^(0)=(1,0; 1,0)ᵀ ergeben die Iterationen unter anderem x^(1)=(0,5000; −0,8636)ᵀ, x^(2)=(0,8494; −0,6413)ᵀ und x^(4)=(0,8127; −0,6646)ᵀ. Ab x^(6) und x^(7) erhält man näherungsweise (0,8122; −0,6650)ᵀ. Der Algorithmus konvergiert damit gegen

x=A⁻¹b ≈ (0,8122; −0,6650)ᵀ.

Konvergenz, Anwendungen und nichtlineare Erweiterung

Das Verfahren konvergiert linear, wenn der Spektralradius der Iterationsmatrix T=−(D+L)⁻¹U kleiner als 1 ist. Der Spektralradius ist der größte Betrag eines Eigenwerts. Je kleiner er ist, desto schneller ist die Konvergenz. Im gegenteiligen Fall divergiert das Verfahren, wenn die rechte Seite einen Anteil eines Eigenvektors zu einem Eigenwert mit Betrag größer als 1 enthält.

Da der Spektralradius praktisch meist schwer zu bestimmen ist, verwendet man hinreichende Kriterien über eine Matrixnorm. Eine besonders gebräuchliche Bedingung ist strikte Diagonaldominanz nach Zeilen:

Σ_(i=1, i≠k)^n |a_ki| < |a_kk| für k=1,…,n.

Die Summe der Beträge aller Nebendiagonalelemente einer Zeile muss also kleiner sein als der Betrag des Diagonalelements. Je größer die kleinste Differenz zwischen rechter und linker Seite dieser Ungleichung ist, desto schneller konvergiert das Verfahren. Zeilen und Spalten können umnummeriert werden, um möglichst große Matrixelemente auf die Diagonale zu bringen; bei Zeilenvertauschungen muss auch die rechte Seite des Gleichungssystems vertauscht werden.

Für große dünnbesetzte Gleichungssysteme aus der Diskretisierung partieller Differentialgleichungen ist das Verfahren ungeeignet. Es wird jedoch als Vorkonditionierer in Krylow-Unterraum-Verfahren oder als Glätter in Mehrgitterverfahren eingesetzt.

Die Idee lässt sich auf nichtlineare Gleichungssysteme f(x)=g übertragen. Im i-ten Schritt wird die i-te Gleichung nach x_i^(k+1) gelöst, während für vorherige Variablen die neuen und für folgende Variablen die bisherigen Werte verwendet werden:

f_i(x_1^(k+1), …, x_(i−1)^(k+1), x_i^(k+1), x_(i+1)^k, …, x_n^k)=g_i.

Dieses Auflösen erfolgt in der Regel selbst durch ein weiteres iteratives Verfahren für nichtlineare Gleichungen. Zur Unterscheidung vom linearen Verfahren heißt dies häufig Gauß-Seidel-Prozess; seine Konvergenz folgt nach dem Banachschen Fixpunktsatz ebenfalls linear.