6.2 Iterative Verfahren
Für eine große Anzahl von Unbekannten ist die Implementierung des Gauß-Algorithmus auf dem Computer nicht geeignet. Einerseits entstehen bedingt durch die hohe Anzahl der Punktoperationen nicht mehr akzeptable Rundungsfehler, andererseits steigt auch der Speicherplatzbedarf an. Der Gauß-Algorithmus hat eine Laufzeitkomplexität von O(n3) für die Erzeugung der Dreiecksmatrix. Für die Rücksubstitution beträgt der Aufwand O(n2).
Iterative Verfahren bieten hier eine Alternative. Der Rechenaufwand ist nur von der Anzahl der Gleichungen und der geforderten Genauigkeit abhängig. Die Laufzeitkomplexität ist linear; sie beträgt k ⋅ O(n), wobei die Konstante k durch die geforderte Rechengenauigkeit festgelegt wird.
Ein besonderer Vorteil von iterativen Algorithmen besteht außerdem darin, dass sie sich relativ einfach implementieren lassen. Allerdings sind iterative Verfahren nicht immer numerisch stabil, d. h., sie konvergieren unter Umständen nicht gegen den Lösungsvektor und liefern dann völlig unbrauchbare Ergebnisse. Um eine numerisch stabile Lösung zu erhalten, müssen die Werte der Hauptdiagonalelemente betragsmäßig größer sein als die übrigen Elemente. Dies ist allerdings nur ein grobes Kriterium; es wird weiter unten präzisiert.
Ich zeige Ihnen nun, wie Sie mit dem Jacobi-Verfahren (Gesamtschrittverfahren) und dem Gauß-Seidel-Verfahren (Einzelschrittverfahren) lineare Gleichungssysteme mit Python lösen können.
6.2.1 Das Jacobi-Verfahren
Die Grundidee ist einfach: Die n-Gleichungen eines linearen Gleichungssystems werden jeweils nach den Unbekannten xi umgestellt. Auf der linken Seite stehen dann die gesuchten Unbekannten. Anschließend wird das System durch Rekursion iterativ gelöst. Als Startwert wird der Nullvektor gewählt.
Betrachten wir das folgende Gleichungssystem mit drei Unbekannten:
Durch die Auflösung nach den Unbekannten ergibt sich die rekursive Gleichung:
Die Iteration wird so lange ausgeführt, bis die gewünschte Genauigkeit erreicht ist. Dies wird über die im Folgenden beschriebene Abbruchbedingung geprüft.
Abbruchbedingungen
Die einfachste Möglichkeit, die Iteration abzubrechen, besteht darin, eine feste Anzahl k = N von Schleifendurchläufen festzulegen.
Alternativ kann man eine bestimmte Genauigkeitsschranke ε vorgeben:
Die Arbeitsweise des Jacobi-Verfahrens soll an dem folgenden Beispiel demonstriert werden:
Die Elemente auf der Hauptdiagonalen sind betragsmäßig größer als die übrigen Elemente. Wenn dieser Fall zutrifft, dann wird die Koeffizientenmatrix als diagonaldominant bezeichnet. Der Begriff diagonaldominant wird in Abschnitt 6.2.3 ausführlich erklärt. Dieses System hat die exakten Lösungen x1 = 1, x2 = 2 und x3 = 3.
Der einfache Lösungsansatz
Listing 6.13 berechnet den Lösungsvektor für das obige Gleichungssystem mit k = 20 Schleifendurchläufen:
01 #13_jacobi1.py
02 import numpy as np
03 N=20
04 a = np.array([[6,2,1],
05 [1,5,3],
06 [2,1,4]],dtype=float)
07 b = np.array([13,20,16],dtype=float)
08 x1=np.empty(N+1,dtype=float)
09 x2=np.empty(N+1,dtype=float)
10 x3=np.empty(N+1,dtype=float)
11 x1[0],x2[0],x3[0]=(0.,0.,0.) #Startwerte
12 print("Jacobi-Verfahren")
13 print("%6s %10s %10s" %("x1","x2","x3"))
14 for k in range(N):
15 x1[k+1]=(b[0]-a[0,1]*x2[k]-a[0,2]*x3[k])/a[0,0]
16 x2[k+1]=(b[1]-a[1,0]*x1[k]-a[1,2]*x3[k])/a[1,1]
17 x3[k+1]=(b[2]-a[2,0]*x1[k]-a[2,1]*x2[k])/a[2,2]
18 print("%2.8f %2.8f %2.8f" %(x1[k],x2[k],x3[k]))
Listing 6.13 Lösung mit dem Jacobi-Verfahren
Ausgabe
x1 x2 x3
0.00000000 0.00000000 0.00000000
2.16666667 4.00000000 4.00000000
0.16666667 1.16666667 1.91666667
1.45833333 2.81666667 3.62500000
0.62361111 1.53333333 2.56666667
1.22777778 2.33527778 3.30486111
0.83743056 1.77152778 2.80229167
1.10910880 2.15113889 3.13840278
0.92655324 1.89513657 2.90766088
1.05034433 2.07009282 3.06293924
0.96614585 1.95216759 2.95730463
1.02306003 2.03238805 3.02888518
0.98438979 1.97805689 2.98037297
1.01058554 2.01489826 3.01329088
0.99281877 1.98990836 2.99098266
1.00486677 2.00684665 3.00611353
0.99669886 1.99535853 2.99585495
1.00223800 2.00314726 3.00281094
0.99848243 1.99786584 2.99809419
1.00102902 2.00144700 3.00129233
Analyse
In den Zeilen 04 bis 06 wird eine diagonaldominante Koeffizientenmatrix a definiert, damit das Verfahren auch konvergiert.
In den Zeilen 08 bis 10 werden drei leere eindimensionale NumPy-Arrays x1, x2 und x3 angelegt. In diesen Arrays sollen die Lösungen des LGS abgespeichert werden.
Innerhalb der for-Schleife (Zeile 14 bis 17) wird der Lösungsvektor rekursiv berechnet. Der jeweils aktuelle Wert wird in Zeile 18 ausgegeben.
An der Ausgabe erkennt man, dass der Lösungsvektor nur langsam konvergiert. Durch die Vergrößerung von N in Zeile 03 können Sie die Genauigkeit der Berechnung noch verbessern.
Übung
Testen Sie das Programm, indem Sie für die Elemente der Hauptdiagonalen kleine bzw. große Zahlen einsetzen, und beobachten Sie das Stabilitätsverhalten. Überprüfen Sie auch, ob die Werte des Inhomogenitätsvektors einen Einfluss auf das Stabilitätsverhalten haben.
Diesen Test sollten Sie mit den drei folgenden Programmvarianten zu den iterativen Verfahren ebenfalls durchführen.
Lösung mit Matrizen
Die erste Programmvariante ist nicht flexibel genug. Wenn mehr als drei Gleichungen gelöst werden sollen, dann müssen jeweils neue Gleichungen in den Schleifenrumpf eingefügt werden. Der Aufwand wäre nicht mehr vertretbar.
Eleganter kann man die iterativen Verfahren implementieren, wenn man die Gleichungssysteme mit Matrizen löst. Dazu betrachten wir ein allgemein formuliertes Gleichungssystem mit drei Unbekannten genauer:
Auf der rechten Seite des Systems ist an den Stellen, wo ursprünglich die Unbekannte xi stand, eine Null eingetragen. Die Faktoren vor den Klammern können als Diagonalmatrix interpretiert werden. In Matrizenschreibweise erhält man dann:
Die Elemente der Hauptdiagonalen bestehen aus den Kehrwerten von aii. Wenn man von den Elementen einer Diagonalmatrix D ihre Kehrwerte berechnet, erhält man die Inverse D–1 dieser Matrix.
Die Koeffizientenmatrix A kann in die strikte linke untere Dreiecksmatrix L und die strikte rechte obere Dreiecksmatrix R zerlegt werden:
Die Vektorschreibweise ergibt dann:
Die Grundidee des Jacobi-Verfahrens beruht also darauf, die Koeffizientenmatrix A des linearen Gleichungssystems in drei Teilmatrizen zu zerlegen: in die Diagonalmatrix D–1, in die strikte linke untere Dreiecksmatrix L und in die strikte rechte obere Dreiecksmatrix R.
Der Term (b – (R + L)) wird in der Fachliteratur auch als Iterationsmatrix B bezeichnet.
Erzeugt wird die strikte linke untere Dreiecksmatrix L mit der NumPy-Funktion tril(a,k=-1). Der Parameter k=-1 ist zwingend notwendig, damit alle Elemente von A oberhalb der Hauptdiagonalen, einschließlich der Hauptdiagonalen selbst, entfernt werden. Das l in tril() steht für lower.
Die strikte rechte obere Dreiecksmatrix R wird mit der NumPy-Funktion triu(a,k=1) erzeugt. Der Parameter k=1 ist auch hier zwingend notwendig, damit alle Elemente von A unterhalb der Hauptdiagonalen, einschließlich der Hauptdiagonalen selbst, entfernt werden. Das u in triu() steht für upper.
Listing 6.14 zeigt die Implementierung des Jacobi-Verfahrens. Damit Sie sich besser vorstellen können, was genau unter einer strikten unteren bzw. strikten oberen Dreiecksmatrix zu verstehen ist, gibt das Programm diese Teilmatrizen mit aus, obwohl diese Ausgabe für die Lösung des Problems nicht notwendig gewesen wäre.
01 #14_jacobi2.py
02 import numpy as np
03 A = np.array([[6,2,1],
04 [1,5,3],
05 [2,1,4]],dtype=float)
06 b = np.array([13,20,16],dtype=float)
07 #Jacobi-Verfahren
08 def jacobi(a,b,N=15):
09 #Dreieckszerlegung
10 L=np.tril(a,k=-1) #links unten
11 R=np.triu(a,k=1) #rechts oben
12 n=np.size(b) #Anzahl der Gleichungen
13 x=np.zeros(n) #Startwerte
14 for k in range(N): #Anzahl der Iterationen
15 for i in range(n):
16 x[i]=(b[i]-(L[i]+R[i])@x)/a[i,i]
17 return x
18 #Ausgabe
19 print("Koeffizientenmatrix\n",A)
20 print("Inhomogenitätsvektor\n",b)
21 print("links-untere-Dreiecksmatrix\n",np.tril(A,k=-1))
22 print("rechts-obere-Dreiecksmatrix\n",np.triu(A,k=1))
23 print("Lösungsvektor\n",jacobi(A,b))
Listing 6.14 Iterative Lösung eines linearen Gleichungssystems mit dem Jacobi-Verfahren
Ausgabe
Koeffizientenmatrix
[[6. 2. 1.]
[1. 5. 3.]
[2. 1. 4.]]
Inhomogenitätsvektor
[13. 20. 16.]
links-untere-Dreiecksmatrix
[[0. 0. 0.]
[1. 0. 0.]
[2. 1. 0.]]
rechts-obere-Dreiecksmatrix
[[0. 2. 1.]
[0. 0. 3.]
[0. 0. 0.]]
Lösungsvektor
[0.99999982 1.99999992 3.00000011]
Analyse
In den Zeilen 14 bis 16 wird das LGS mit dem Vektoransatz des Jacobi-Verfahrens gelöst. Schon nach 15 Schleifendurchläufen der äußeren Schleife erzielt das Programm eine akzeptable Genauigkeit. Zeile 16 zeigt deutlich, dass keine Inverse berechnet werden muss, denn der rechtsseitige Term (b[i]-(L[i]+R[i])@x) wird nur durch die Diagonalelemente a[i,i] elementweise dividiert.
6.2.2 Das Gauß-Seidel-Verfahren
Dieses Verfahren wurde zwar von Gauß schon 1823 entwickelt, aber erst im Jahre 1874 von Ludwig Seidel (1821–1896) veröffentlicht.
Die Grundidee ist einfach. Man löst wie beim Jacobi-Verfahren die einzelnen Gleichungen eines linearen Gleichungssystems jeweils nach den Unbekannten x1, x2 … xn auf, legt Startwerte fest und berechnet dann schrittweise die Lösungen. Im Unterschied zum Jacobi-Verfahren werden jedoch die in einem Schleifendurchlauf zuvor berechneten Werte in die rechte Seite des rekursiven Gleichungssystems eingesetzt. Es wird also mit bereits bekannten Werten weitergerechnet.
Der einfache Lösungsansatz
Der Lösungsansatz gleicht dem Jacobi-Verfahren. Der Unterschied besteht darin, bereits berechnete Werte xk+1 in die rechte Seite der Gleichung einzusetzen. Ein einfacher Lösungsansatz lässt sich wie folgt formulieren:
Listing 6.15 zeigt, wie das Gleichungssystem
mit dem Gauß-Seidel-Verfahren gelöst werden kann. Eigentlich wäre es überflüssig gewesen, für die drei Unbekannten x1, x2 und x3 jeweils ein Array anzulegen (Zeilen 08 bis 10), wie es beim Jacobi-Verfahren noch zwingend erforderlich war. Damit aber der Unterschied zwischen den beiden Verfahren deutlich wird, habe ich diesen Weg gewählt.
Für die Wiederholstruktur (Zeilen 15 bis 17) könnte eine einfachere Variante wie folgt implementiert werden:
x1=(b[0]-a[0,1]*x2-a[0,2]*x3)/a[0,0]
x2=(b[1]-a[1,0]*x1-a[1,2]*x3)/a[1,1]
x3=(b[2]-a[2,0]*x1-a[2,1]*x2)/a[2,2]
Das Programm berechnet mit nur 15 Iterationen bereits akzeptable Ergebnisse.
01 #15_seidel1.py
02 import numpy as np
03 N=15
04 a = np.array([[6,2,1],
05 [1,5,3],
06 [2,1,4]],dtype=float)
07 b = np.array([13,20,16],dtype=float)
08 x1=np.empty(N+1,dtype=float)
09 x2=np.empty(N+1,dtype=float)
10 x3=np.empty(N+1,dtype=float)
11 x1[0],x2[0],x3[0]=(0.,0.,0.) #Startwerte
12 print("Gauss-Seidel-Verfahren")
13 print("%6s %10s %10s" %("x1","x2","x3"))
14 for k in range(N):
15 x1[k+1]=(b[0]-a[0,1]*x2[k]-a[0,2]*x3[k])/a[0,0]
16 x2[k+1]=(b[1]-a[1,0]*x1[k+1]-a[1,2]*x3[k])/a[1,1]
17 x3[k+1]=(b[2]-a[2,0]*x1[k+1]-a[2,1]*x2[k+1])/a[2,2]
18 print("%2.8f %2.8f %2.8f" %(x1[k+1],x2[k+1],x3[k+1]))
Listing 6.15 Lösung eines linearen Gleichungssystems mit dem Gauß-Seidel-Algorithmus
Ausgabe
Gauss-Seidel-Verfahren
x1 x2 x3
2.16666667 3.56666667 2.02500000
0.64027778 2.65694444 3.01562500
0.77841435 2.03494213 3.10205729
0.97134307 1.94449701 3.02820421
1.01380029 1.98031741 2.99802050
1.00689078 1.99980955 2.99660222
1.00062978 2.00191271 2.99920693
0.99949461 2.00057692 3.00010847
0.99978962 1.99997700 3.00011094
0.99998918 1.99993560 3.00002151
1.00001788 1.99998352 2.99999518
1.00000630 2.00000163 2.99999644
1.00000005 2.00000212 2.99999944
0.99999938 2.00000046 3.00000019
0.99999982 1.99999992 3.00000011
Analyse
Der Algorithmus hat eine ähnliche Struktur wie das Jacobi-Verfahren, jedoch mit dem gravierenden Unterschied, dass unterhalb der Hauptdiagonalen die schon bereits berechneten Linkswerte x[k+1] stehen. Das Gauß-Seidel-Verfahren konvergiert deshalb schneller als das Jacobi-Verfahren. Schon nach 15 Iterationen erzielt das Programm eine zufriedenstellende Genauigkeit.
Gauß-Seidel-Verfahren mit Matrizen
Betrachten wir den einfachen Lösungsansatz genauer, indem wir auf der Hauptdiagonalen Nullen eintragen:
Jetzt erkennen wir, dass unterhalb der Hauptdiagonalen die schon berechneten Werte mit den Indexen k + 1 stehen. Die Koeffizientenmatrix wird wieder in die strikte untere linke Matrix L und die strikte obere rechte Matrix R zerlegt:
In Vektorschreibweise ergibt sich dann:
Alle Terme mit dem Index k + 1 müssen auf die linke Seite der Gleichung gebracht werden:
Durch Auflösen nach xk+1 ergibt sich dann:
Mit der linken unteren Dreiecksmatrix
erhält man den iterativen Algorithmus in Vektorschreibweise
für das Gauß-Seidel-Verfahren.
Diese Matrizengleichung lässt sich besonders einfach programmieren. Sie benötigt nur eine Programmzeile. Die linke untere Dreiecksmatrix wird mit der NumPy-Funktion tril(A) erzeugt.
Listing 6.16 zeigt die Umsetzung:
01 #16_seidel2.py
02 import numpy as np
03 from numpy.linalg import inv
04 A = np.array([[6,2,1],
05 [1,5,3],
06 [2,1,4]],dtype=float)
07 b = np.array([13,20,16],dtype=float)
08 #Gauss-Seidel-Verfahren
09 def seidel(A, b, N=15):
10 L=np.tril(A)
11 invL=inv(L)
12 R=A-L
13 n=np.size(b)
14 x=np.zeros(n) #Startwerte
15 for k in range(N):
16 x = invL@(b - R@x)
17 #print(x)
18 return x
19 #Ausgabe
20 print("Lösungsvektor\n",seidel(A,b))
Listing 6.16 Algorithmus für das Gauß-Seidel-Verfahren
Ausgabe
Lösungsvektor
[0.99999982 1.99999992 3.00000011]
Analyse
Die eigentliche Berechnung des Lösungsvektors x benötigt nur eine Programmzeile (Zeile 16). Beachten Sie allerdings, dass die Inverse der unteren linken Dreiecksmatrix berechnet werden muss (Zeile 11).
Wenn Sie das Konvergenzverhalten beobachten wollen, dann können Sie den Kommentar in Zeile 17 entfernen.
6.2.3 Konvergenzverhalten
Beide Verfahren konvergieren linear. Das Gauß-Seidel-Verfahren konvergiert jedoch schneller als das Jacobi-Verfahren. Sie konvergieren aber nur, wenn die Beträge der Diagonalelemente größer sind als die Zeilensummen der Koeffizientenmatrix. Man spricht in diesem Fall davon, dass die Matrix diagonaldominant ist. Neben diesen sogenannten Zeilensummenkriterien gibt es noch das Spaltensummenkriterium, das hier aber nicht behandelt werden soll.
Wenn eine Koeffizientenmatrix nicht diagonaldominant ist, dann liefern die hier vorgestellten iterativen Verfahren völlig falsche Ergebnisse.
Das Zeilensummenkriterium
ist eine hinreichende Bedingung für die numerische Stabilität des Jacobi- und des Gauß-Seidel-Verfahrens. Wenn die Koeffizientenmatrix A diagonaldominant ist, dann konvergiert sowohl das Jacobi- als auch das Gauß-Seidel-Verfahren für Ax = b [Knorrenschild: 65].
Tabelle 6.2 zeigt ein Beispiel für numerisch instabiles Verhalten. Diese Matrix ist nicht diagonaldominant.
|
Zeile |
Koeffizientenmatrix |
Zeilensummen |
Dialogelemente |
||||
|---|---|---|---|---|---|---|---|
|
1 |
3 |
2 |
1 |
3 |
= |
3 |
|
|
2 |
1 |
2 |
3 |
4 |
> |
2 |
|
|
3 |
2 |
1 |
4 |
3 |
< |
4 |
|
Tabelle 6.2 Zeilensummenkriterium für ein numerisch instabiles Verhalten
Tabelle 6.3 zeigt ein Beispiel für numerisch stabiles Verhalten.
|
Koeffizientenmatrix |
Zeilensummen |
Dialogelemente |
|||||
|---|---|---|---|---|---|---|---|
|
1 |
6 |
2 |
1 |
3 |
< |
6 |
|
|
2 |
1 |
5 |
3 |
4 |
< |
5 |
|
|
3 |
2 |
1 |
4 |
3 |
< |
4 |
|
Tabelle 6.3 Zeilensummenkriterium für ein numerisch stabiles Verhalten
Auf den ersten Blick fällt sofort auf, dass die Diagonalelemente betragsmäßig größer sind als die Summe der übrigen Elemente ihrer jeweiligen Zeile. Diese Matrix ist also diagonaldominant.
Die Diagonaldominanz ist nur ein hinreichendes Kriterium für die numerische Stabilität des Jacobi- und Gauß-Seidel-Verfahrens. Es gibt auch Ausnahmen, bei denen beide Verfahren mit nicht diagonaldominanten Koeffizientenmatrizen konvergieren. So konvergiert z. B. das Gauß-Seidel-Verfahren für alle symmetrischen positiv definiten Matrizen A [Knorrenschild: 64].
Eine Matrix A ist positiv definit, wenn gilt [Strang: 348]:
Ein notwendiges und hinreichendes Konvergenzkriterium liefert der Spektralradius ρ der Iterationsmatrix B. Er muss kleiner 1 sein: ρ(B) < 1
Spektralradius
Der Spektralradius ρ einer quadratischen Matrix B ist der Betrag des betragsmäßig größten Eigenwerts |λ|max von B. Das heißt, er ist definiert durch:
Wie Eigenwerte berechnet werden, beschreibt [Strang: 289] ausführlich.
Das Jacobi-Verfahren hat die Iterationsmatrix
und für das Gauß-Seidel-Verfahren gilt:
Obwohl das Gauß-Seidel-Verfahren wegen des besseren Konvergenzverhaltens dem Jacobi-Verfahren überlegen zu sein scheint, gibt es Linearsysteme, für die das Jacobi-Verfahren konvergiert, das Gauß-Seidel-Verfahren jedoch nicht. Umgekehrt gibt es Fälle, bei denen nur das letztere Verfahren konvergiert. Wenn die Koeffizientenmatrix A diagonaldominant ist, dann konvergieren beide Verfahren für einen beliebigen Inhomogenitätsvektor b [Faires: 310].