Lineare Gleichungssysteme

Lineare Gleichungssysteme treten oft als Teilprobleme bei numerischen Verfahren auf, z.B.:

• Lineare Ausgleichsprobleme (Normalgleichungen)
• Partielle Differentialgleichungen

  • bei manchen impliziten Verfahren muß in jedem Zeitschritt ein
    lineares Gleichungssystem gelöst werden
  • Randwertprobleme bei elliptischen Systemen (z.B. Potentialgleichung)
  • man hat oft spezielle, dünnbesetzte (sparse ) Matrizen (Bandmatrizen )

Direkte Verfahren
Lösen das Problem durch geeignete Zerlegung der Koeffizientenmatrix nach endlich vielen Schritten von algebraischen Operationen, sodaß kein Verfahrensfehler auftritt. Abgesehen von Rundungsfehlern erhält man die exakte Lösung. Allerdings können wegen der großen Anzahl von Punktoperationen µ n³ die akkumulierten Rundungsfehler das Ergebnis erheblich verfälschen, sodaß die Lösung völlig unbrauchbar werden kann.
Beispiele: Gaußscher Algorithmus, LU-Zerlegung.

Iterationsverfahren
Ausgehend von einer Anfangsnäherung für die Lösung (Startvektor) wird diese schrittweise durch Iteration verbessert. Obwohl hier Verfahrens- und Rundungsfehler auftreten, sind Iterationsverfahren vor allem für große Systeme mit schwach besetzten Matrizen geeignet (pro Iterationsschritt ist die Anzahl der Punktoperationen µ n²).
Beispiele: Jacobi, Gauß-Seidel, SOR, ADI, CG.

Das Gaußsche Eliminationsverfahren

Wir betrachten ein System von n linearen Gleichungen in n Unbekannten x₁, ¼, x_n :

a11 x1 + a12 x2 +  ¼  + a1n xn = b1 a21 x1 + a22 x2 +  ¼  + a2n xn = b2 : : : : an1 x1 + an2 x2 +  ¼  + ann xn = bn     ü ï ï ý ï ï þ

(1)

oder in Matrixform:

A = æ ç ç ç è a11 ¼ a1n : : an1 ¼ ann ö ÷ ÷ ÷ ø   ,

x = æ ç ç ç è x1 : xn ö ÷ ÷ ÷ ø   ,

b = æ ç ç ç è b1 : bn ö ÷ ÷ ÷ ø   ,

A · x = b

Beispiel zum Gaußschen Eliminationsverfahren

2 x1 + x2 + x3 = 5        G1(0) 4 x1 - 6 x2 = -2        G2(0) -2 x1 + 7 x2 + 2 x3 = 9        G3(0) 2 x1 + x2 + x3 = 5        G1(0) - 8 x2 - 2 x3 = -12        G2(1) = G2(0) - 2 ·G1(0) 8 x2 + 3 x3 = 14        G3(1) = G3(0) - ( -  G1(0)) 2 x1 + x2 + x3 = 5        G1(0) - 8 x2 - 2 x3 = -12        G2(1) x3 = 2        G3(2) = G3(1) - ( -  G2(1))

Damit ist das ursprüngliche Gleichungssystem in ein äquivalentes Gleichungssystem von oberer Dreiecksgestalt übergeführt, welches sich von rückwärts her auflösen läßt:

Aus der letzten Gleichung erhält man

x3 = 2     .

Einsetzen von x₃ in die vorletzte Gleichung liefert

x2 = (-12 + 4)/(-8) = 1     .

Aus Gleichung G1⁽⁰⁾ erhält man schließlich

x1 = (5 - 1 - 2)/2 = 1     .

Die Lösung x der ursprünglichen Gleichung lautet also

x   =   æ ç ç ç è x1 x2 x3 ö ÷ ÷ ÷ ø   =   æ ç ç ç è 1 1 2 ö ÷ ÷ ÷ ø     .

Die Koeffizientenmatrix A sei nicht singulär, det A ¹ 0. Dann ist (1) eindeutig lösbar. Das Prinzip des Gaußschen Algorithmus besteht nun darin, (1) in ein Gleichungssystem von oberer Dreiecksgestalt (mit einer oberen Dreiecksmatrix U als neuer Koeffizientenmatrix) überzuführen, welches dieselbe Lösung wie (1) besitzt:

u11 x1 + u12 x2 + u13 x3 + ¼ + u1n xn = b¢1 u22 x2 + u23 x3 + ¼ + u2n xn = b¢2 ··· : : un-1,n-1 xn-1 + un-1,n xn = b¢n-1 unn xn = b¢n    ü ï ï ï ý ï ï ï þ

(2)

Falls alle u_ii ¹ 0 sind, läßt sich (2), angefangen von der letzten Gleichung, von rückwärts her sukzessive auflösen:

xn = b¢n  / unn xn-1 = ( b¢n-1 - un-1,n xn )  / un-1,n-1 : x2 = ( b¢2 - u23 x3 - ¼- u2n xn )  / u22 x1 = ( b¢1 - u12 x2 - u13 x3 - ¼- u1n xn )  / u11

Allgemein erhält man für das Rückwärtseinsetzen (i = n-1, n-2, ¼, 1):

xn = b¢n  / unn xi = (  b¢i  -  n å j = i + 1  uij  xj  )  / uii

(3)

Alle nicht singulären Gleichungssystme lassen sich auf die Form (2) bringen: Die Lösungsmenge eines linearen Gleichungssystems ändert sich nämlich nicht, wenn man

• die Reihenfolge der Gleichungen vertauscht
• zu einer Gleichung das Vielfache einer anderen Gleichung addiert
  und eine der beiden Gleichungen durch diese Summe ersetzt.

Im ersten Schritt definieren wir n - 1 Multiplikatoren für die Zeilen k = 2 bis k = n,

m(1)k  =   a(0)k1 a(0)11        k = 2, 3, ¼, n     ,     a(0)11 ¹ 0     ,

multiplizieren die erste Gleichung mit m⁽¹⁾₂ und subtrahieren sie von der zweiten, dann multiplizieren wir die erste Gleichung mit m⁽¹⁾₃ und subtrahieren sie von der dritten usw. Dadurch wird x₁ von der 2. bis zur n-ten Gleichung eliminiert:

a(0)11 x1 + a(0)12 x2 + ¼ + a(0)1n xn = b(0)1 0 a(1)22 x2 + ¼ + a(1)2n xn = b(1)2 : : : : 0 a(1)n2 x2 + ¼ + a(1)nn xn = b(1)n        ü ï ï ï ý ï ï ï þ

(4)

Die erste Zeile bleibt dabei unverändert und heißt Pivotzeile . Das aktuell zur Elimination benutzte Diagonalelement a⁽⁰⁾₁₁ heißt Pivotelement (frz. pivot: Drehpunkt, Angelpunkt). Die neuen Koeffizienten a⁽¹⁾_kj (j = 2, ¼, n) und rechten Seiten b⁽¹⁾_k der Zeilen k = 2, 3, ¼, n lauten

a(1)kj = a(0)kj  - m(1)k a(0)1j b(1)k = b(0)k  - m(1)k  b(0)1

Dieses Verfahren wird nun mit den verbleibenden (n - 1) Gleichungen wiederholt, indem man wieder Multiplikatoren

m(2)k  =   a(1)k2 a(1)22        k = 3, 4, ¼, n     ,     a(1)22 ¹ 0     ,

definiert, die zweite Gleichung des Systems (4) als neue Pivotzeile verwendet, mit dem entsprechenden m⁽²⁾_k multipliziert, und sie dann von der k-ten Gleichung subtrahiert.

Im i-ten Schritt definiert man die Multiplikatoren

m(i)k  =   a(i-1)ki a(i-1)ii              k = i + 1, ¼, n     ,     a(i-1)ii ¹ 0     ,

(5)

und subtrahiert von der k-ten Gleichung (k = i + 1, ¼, n) das m⁽ⁱ⁾_k-fache der i-ten Gleichung (der aktuellen Pivotzeile). Die neuen Koeffizienten a⁽ⁱ⁾_kj (j = i + 1, ¼, n) und rechten Seiten b⁽ⁱ⁾_k der Zeilen k = i + 1, ¼, n sind dann

a(i)kj = a(i-1)kj  - m(i)k a(i-1)ij b(i)k = b(i-1)k  - m(i)k  b(i-1)i

(6)

Nach dem i-ten Schritt hat das System die Gestalt

a(0)11 x1 + a(0)12 x2 + a(0)13 x3 + a(0)14 x4 + ¼ + a(0)1n xn = b(0)1 a(1)22 x2 + a(1)23 x3 + a(1)24 x4 + ¼ + a(1)2n xn = b(1)2 ··· : : a(i-1)ii xi + a(i-1)i,i+1 xi+1 + ¼ + a(i-1)in xn = b(i-1)i 0 a(i)i+1,i+1 xi+1 + ¼ + a(i)i+1,n xn = b(i)i+1 : : : : 0 a(i)n,i+1 xi+1 + ¼ + a(i)nn xn = b(i)n

Damit wird x_i aus allen Zeilen k = i + 1, ¼, n eliminiert und die ersten i Spalten enthalten unterhalb der Hauptdiagonale nur mehr Nullen.
Nach n - 1 Schritten hat man schließlich die Form (2) mit einer oberen Dreiecksmatrix als Koeffizientenmatrix erreicht:

a(0)11 x1 + a(0)12 x2 + ¼ + a(0)1n xn = b(0)1 a(1)22 x2 + ¼ + a(1)2n xn = b(1)2 ··· : : a(n-2)n-1,n-1 xn-1 + a(n-2)n-1,n xn = b(n-2)n-1 a(n-1)nn xn = b(n-1)n    ü ï ï ï ý ï ï ï þ

(7)

Die oberen Indizes (i) geben an, wie oft die Koeffizienten a⁽ⁱ⁾_kj bzw. Elemente b⁽ⁱ⁾_k der rechten Seite während der Eliminationsrechnung höchstens verändert wurden. Die Gesamtanzahl der Multiplikationen und Divisionen bei der Elimination und beim Rückwärtseinsetzen ist ( n³ + 3 n² - n ) / 3  = O(n³). Der Zeitaufwand für die Lösung eines allgemeinen linearen Gleichungssystems steigt beim Gaußverfahren also mit der dritten Potenz der Anzahl n der Gleichungen: Bei einer Verdopplung von n erhöht sich die Rechenzeit auf ca. das Achtfache. Wegen der großen Anzahl von Rechenoperationen besteht die Gefahr, daß sich Rundungsfehler akkumulieren. Man sollte daher die Rechnungen nur in doppelter Genauigkeit (double ) durchführen.

Um unnötigen Rechenaufwand (und damit eine Anhäufung von Rundungsfehlern) zu vermeiden, sollte man auf keinen Fall ein Gleichungssystem A ·x = b lösen, indem man zuerst die zu A inverse Matrix A^-1 berechnet (dies erfordert die Lösung von n linearen Gleichungssystemen mit den n Spalten der (n ×n)-Einheitsmatrix als rechte Seiten) und dann durch Multiplikation die Lösung x = A^-1 ·b ermittelt!

Durchführbarkeit des Gaußschen Eliminationsverfahrens

Bei den Eliminationsschritten wurde vorausgesetzt, daß die jeweiligen Diagonalelemente (Pivotelemente) a^(i-1)_ii ¹ 0 sind. Ist das nicht der Fall, wie zum Beispiel bei

x2 = 1 x1 + x2 = 2   ,

so versagt der Gaußsche Algorithmus in seiner bisherigen Form. Wenn das System (1) nicht singulär ist, kann man durch Zeilenvertauschungen immer erreichen, daß das gerade aktuelle Pivotelement a^(i-1)_ii ¹ 0 wird. Man wählt jedoch als neue Pivotzeile im i-ten Schritt nicht die erstbeste Zeile des verbleibenden Systems mit a^(i-1)_ki ¹ 0 (k = i, ¼, n), sondern diejenige Zeile k₀ mit dem betragsgrößten Element in Spalte i (Abb. 1):

| a(i-1)k0i |  =  max k = i, ¼, n  | a(i-1)ki |

Die so gefundene neue Pivotzeile k₀ wird dann mit der aktuellen Pivotzeile i vertauscht. Damit bleiben die Multiplikatoren |m⁽ⁱ⁾_k| £ 1 und Rundungsfehler werden besonders klein gehalten. Dieses Verfahren heißt Spaltenpivotsuche oder teilweise Pivotsuche :

• Die Zeilenvertauschungen brauchen in einem Programm nicht tatsächlich ausgeführt zu werden, es genügt eine Indexvertauschung: Man definiert einen zusätzlichen Merkvektor (Permutationsvektor ) p[1],..,p[n], in den man alle Vertauschungen einträgt. Er wird zu Beginn mit p[1]=1,..,p[n]=n initialisiert und jeder Zeilenindex i wird durch p[i] ersetzt, z.B: a[i][j] ® a[p[i]][j].
• Ist A diagonaldominant, |a_ii| > å_{j=1, j ¹ i}ⁿ |a_ij| (i = 1, ¼, n), so ist das Gaußverfahren ohne Pivotsuche durchführbar. Um Rundungsfehler zu minimieren, sollte aber eine Pivotsuche auch dann durchgeführt werden, wenn keines der aktuellen Pivotelemente Null ist!
• In der Praxis kann es vorkommen, daß sich die Koeffizienten a_ij um Größenordnungen unterscheiden. Damit in solchen Fällen die Pivotsuche trotzdem wirkungsvoll ist, führt man eine sogenannte skalierte Spaltenpivotsuche durch: Man definiert für jede Zeile k einen Skalierungsfaktor s_k : = max_j=1,¼,n | a_kj | und verwendet zur Bestimmung der Pivotzeile im i-ten Schritt max_k=i,¼,n ( | a^(i-1)_ki | / s_k ). Damit wird sichergestellt, daß das betragsgrößte Element jeder Zeile die relative Größe von 1 besitzt. Die Skalierung wird in einem Programm nur zum Vergleich bei der Suche nach Pivotzeilen durchgeführt; die Koeffizienten und rechten Seiten selbst bleiben dabei unverändert, um keine zusätzlichen Rundungsfehler zu erzeugen. Zeilenvertauschungen müssen im Vektor s[] der Skalierungsfaktoren s_k mitberücksichtigt werden: s[i] ® s[p[i]].

Ein singuläres Gleichungssystem wird vom Gaußschen Algorithmus daran erkannt, daß an einer bestimmten Stelle die Pivotsuche erfolglos ist, z.B. ist dann im i-ten Schritt a^(i-1)_ki = 0 für alle i £ k £ n. Durch Rundungsfehler kann es aber passieren, daß ein singuläres System als lösbar erscheint. In der Praxis erklärt man daher ein lineares Gleichungssystem für nicht behandelbar, wenn ein Pivotelement betragsmäßig kleiner als ein vorgegebenes e wird (z.B. e = 10^-8). Die Matrix A ist dann numerisch singulär .

Abb. 1:  Spaltenpivotsuche im i-ten Schritt.

Abb. 2:  Gut und schlecht konditioniertes 2×2-System.

LU-Zerlegung

Bei diesem Verfahren zur Lösung von (1) wird die Koeffizientenmatrix A zunächst in ein Produkt von zwei Dreiecksmatrizen zerlegt,

A  = L ·U     ,

wobei L eine untere (lower ) Dreiecksmatrix und U eine obere (upper) Dreiecksmatrix ist:

L  =  æ ç ç ç ç ç è l11 0 0 ¼ 0 l21 l22 0 ¼ 0 : : ··· : ln1 ln2 ln3 ¼ lnn ö ÷ ÷ ÷ ÷ ÷ ø   ,              U  =  æ ç ç ç ç ç è u11 u12 u13 ¼ u1n 0 u22 u23 ¼ u2n : : ··· : 0 0 0 ¼ unn ö ÷ ÷ ÷ ÷ ÷ ø     .

Schreibt man nun das ursprüngliche Gleichungssystem A ·x  = b als

L ·( U ·x )  = b   ,

so gewinnt man die Lösung von (1) in einem zweistufigen Vorgang: Man löst zuerst

L ·y  = b

(8)

nach y und verwendet den Lösungsvektor y als rechte Seite in

U ·x  = y   .

(9)

Wegen der Dreiecksform von L und U sind die beiden letzten Gleichungen leicht zu lösen. Man bestimmt zunächst den Hilfsvektor y in (8) durch Vorwärtseinsetzen (i = 2, ¼, n):

y1 = b1  / l11 yi = (  bi  -  i - 1 å j = 1  lij  yj  )  / lii

(10)

Den gesuchten Lösungsvektor x erhält man dann aus (9) durch Rückwärtseinsetzen wie im Gaußverfahren, vgl. (3). Vorteile der LU-Aufspaltung:

• Kennt man eine Faktorisierung A = L ·U, so läßt sich die Lösung von A ·x = b in O(n²) Rechenoperationen gewinnen.
• Die Faktoren L und U sind unabhängig von b und können, einmal bestimmt, für verschiedene rechte Seiten b verwendet werden.
• Als Nebenprodukt einer LU-Zerlegung erhält man det A = (-1)^k u₁₁ u₂₂ ¼u_nn, wobei k die Anzahl der Zeilenvertauschungen ist (siehe unten).

Wenn das Gaußsche Eliminationsverfahren ohne Pivotsuche durchführbar ist, dann liefert es eine Dreieckszerlegung von A,

A = L ·U   ,

mit

L  =  æ ç ç ç ç ç ç ç è 1 0 ¼ 0 0 l21 1 ¼ 0 0 : ··· : ln-1,1 ln-1,2 ¼ 1 0 ln1 ln2 ¼ ln,n-1 1 ö ÷ ÷ ÷ ÷ ÷ ÷ ÷ ø     und     U  =  æ ç ç ç ç ç è a(0)11 a(0)12 ¼ a(0)1n 0 a(1)22 ¼ a(1)2n : ··· : 0 0 ¼ a(n-1)nn ö ÷ ÷ ÷ ÷ ÷ ø   ,

wobei die l_ki genau die im i-ten Schritt definierten Multiplikatoren m⁽ⁱ⁾_k sind, vgl. (5), und U die nach (n - 1) Eliminationsschritten erhaltene Koeffizientenmatrix A^(n-1) ist, vgl. (7). Die in L auftretenden Zahlen l_ki unterhalb der Hauptdiagonale können in dem (für die Elimination nicht benötigten) unteren Teil von A⁽ⁱ⁾ gespeichert werden (die l_ii = 1 werden nicht abgespeichert); damit wird die untere Dreiecksmatrix L im Laufe der Elimination spaltenweise aufgebaut.

Ist das Gaußsche Eliminationsverfahren nur mit Pivotsuche durchführbar, so liefert es eine Zerlegung von P ·A,

P ·A  = L ·U   ,

wo P die Permutationsmatrix ist, welche die im Laufe des Verfahrens durchgeführten Zeilenvertauschungen beschreibt.

Beispiel: Multiplikation einer (3×3)-Matrix A mit

P  =  æ ç ç ç è 0 1 0 1 0 0 0 0 1 ö ÷ ÷ ÷ ø

vertauscht die erste und zweite Zeile von A.

Man hat also in diesem Fall ( P ·A ) ·x  = P ·b  =: c und löst

L ·U ·x  = c

wie oben mit Vorwärts- und Rückwärtseinsetzen. Programmiertechnisch bedeutet das, daß man dabei auf die Zeilenindizes der ermittelten LU-Matrix und der originalen rechten Seite über den Permutationsvektor p[] zugreift; dazu muß während der Elimination die untere Dreiecksmatrix L in der ursprünglichen Zeilenordnung gespeichert werden, also a[p[k]][j] = zmult, und beim Vorwärts- und Rückwärtseinsetzen muß p[] zur Verfügung stehen.

Neben dem Gaußschen Algorithmus gibt es noch direkte Verfahren zur Berechnung von L und U (Crout, Doolittle, Banachiewicz), welche nur halb so viele Rechenoperationen benötigen. Allerdings ist hier die Pivotierung komplizierter als beim Gaußverfahren.

Fehler und Kondition

Die mit Hilfe direkter Methoden ermittelte Lösung x^* eines linearen Gleichungssystems ist meist nicht die exakte Lösung x , da (1) Rundungsfehler akkumuliert werden, und da (2) Ungenauigkeiten in den Ausgangsgrößen A und b auch Ungenauigkeiten in der Lösung hervorrufen können. Leider ist das Residuum r   :=  b - A ·x^* als Maß für die Güte einer Näherungslösung x^* wenig geeignet: Ist nämlich die Länge || r || des Residuums klein, so kann man daraus nicht schließen, daß auch die Länge || x - x^* || des unbekannten Fehlers x - x^* klein ist. Zur Beurteilung der Güte von x^* ist also eine Probe, d.h. das Einsetzen der Näherungslösung x^* in das gegebene Gleichungssystem A ·x = b, nicht hinreichend.

Beispiel: Das System

æ ç è 0 1 a -1 ö ÷ ø · æ ç è x y ö ÷ ø    =    æ ç è 0 0 ö ÷ ø

beschreibt den Schnittpunkt der Geraden y = 0 mit der Geraden y = ax und hat die exakte Lösung x = (0, 0)^T. Für ein beliebiges x^* = (x^*, 0)^T ist das Residuum r = ( 0, - ax^* )^T. Ist etwa x^* = 10⁵ , a = 10^-10 , so ist zwar ||r|| = ax^* » 10^-5, trotzdem ist der Fehler || x - x^* || = 10⁵. Die Lösung x^* hat also mit der exakten Lösung fast nichts zu tun, produziert aber beim Einsetzen nur ein kleines Residuum.

Es stellt sich heraus, daß die sogenannte Kondition k(A) : = ||A||  ||A^-1 || ³ 1 der Koeffizientenmatrix A entscheidend dafür ist, ob man bei der numerischen Lösung des linearen Gleichungssystems mit Schwierigkeiten zu rechnen hat. k(A) ist der Faktor, um den die relativen Fehler in A und b zum relativen Fehler der Lösung verstärkt werden:

Wenn kleine Änderungen in den Ausgangsdaten A und b große Änderungen in der Lösung x^* hervorrufen, wenn also k(A) groß ist, so spricht man von einem schlecht konditionierten System.

Das Paradebeispiel für eine schlecht konditionierte Matrix ist die sogenannte Hilbertmatrix

Hn   :=   æ ç ç ç ç ç ç ç ç ç ç ç è 1  1 2 ¼  1 n  1 2  1 3 ¼  1 n + 1 : :  1 n  1 n + 1 ¼  1 2n - 1 ö ÷ ÷ ÷ ÷ ÷ ÷ ÷ ÷ ÷ ÷ ÷ ø    =   ( hij )    =    æ è  1 i + j - 1 ö ø     für     i,j = 1, ¼, n   .

Obwohl diese Matrix zunächst gutartig aussieht (symmetrisch, positiv definit, regulär), sind Gleichungssysteme mit dieser Matrix nur schwer zu lösen: k(H_n) wächst exponentiell mit n. Zum Beispiel hat das System H_n ·x = b für b_i = åⁿ_j=1 h_ij die exakte Lösung x = (1, ¼, 1)^T. Ein durch Rundungsfehler leicht verändertes System hat jedoch eine Lösung, die weit von x abweicht.