Lineare Gleichungssysteme

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


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 n3 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 n2).
Beispiele: Jacobi, Gauß-Seidel, SOR, ADI, CG.



Das Gaußsche Eliminationsverfahren

Wir betrachten ein System von n linearen Gleichungen in n Unbekannten x1, , xn :
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 x3 in die vorletzte Gleichung liefert
x2 = (-12 + 4)/(-8) = 1     .
Aus Gleichung G1(0) 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
=
b1
u22 x2
+
u23 x3
+
+
u2n xn
=
b2
···
:
:
un-1,n-1 xn-1
+
un-1,n xn
=
bn-1
unn xn
=
bn
  







(2)
Falls alle uii 0 sind, läßt sich (2), angefangen von der letzten Gleichung, von rückwärts her sukzessive auflösen:
xn
=
bn  / unn
xn-1
=
( bn-1 - un-1,n xn )  / un-1,n-1
:
x2
=
( b2 - u23 x3 - - u2n xn )  / u22
x1
=
( b1 - u12 x2 - u13 x3 - - u1n xn )  / u11

Allgemein erhält man für das Rückwärtseinsetzen (i = n-1, n-2, , 1):
xn
=
bn  / unn
xi
=
(  bi  -  n

j = i + 1 
uij  xj  )  / uii
(3)

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

Wir verwenden nun diese Eigenschaften, um das ursprüngliche Gleichungssystem (1) so umzuformen, daß alle Elemente von A unterhalb der Hauptdiagonale verschwinden. In der Ausgangsform bezeichnen wir das gegebene Gleichungssystem (1) mit
A(0) ·x  = b(0)   .

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(1)2 und subtrahieren sie von der zweiten, dann multiplizieren wir die erste Gleichung mit m(1)3 und subtrahieren sie von der dritten usw. Dadurch wird x1 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(0)11 heißt Pivotelement (frz. pivot: Drehpunkt, Angelpunkt). Die neuen Koeffizienten a(1)kj (j = 2, , n) und rechten Seiten b(1)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(2)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(i)k-fache der i-ten Gleichung (der aktuellen Pivotzeile). Die neuen Koeffizienten a(i)kj (j = i + 1, , n) und rechten Seiten b(i)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 xi 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(i)kj bzw. Elemente b(i)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 ( n3 + 3 n2 - n ) / 3  = O(n3). 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 k0 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 k0 wird dann mit der aktuellen Pivotzeile i vertauscht. Damit bleiben die Multiplikatoren |m(i)k| 1 und Rundungsfehler werden besonders klein gehalten. Dieses Verfahren heißt Spaltenpivotsuche oder teilweise Pivotsuche :


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 .


pivot2.png
Abb. 1:  Spaltenpivotsuche im i-ten Schritt.



cond1.png
        cond2.png
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:

Wie findet man nun die Matrizen L und U ?

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 lki genau die im i-ten Schritt definierten Multiplikatoren m(i)k sind, vgl. (5), und U die nach (n - 1) Eliminationsschritten erhaltene Koeffizientenmatrix A(n-1) ist, vgl. (7). Die in L auftretenden Zahlen lki unterhalb der Hauptdiagonale können in dem (für die Elimination nicht benötigten) unteren Teil von A(i) gespeichert werden (die lii = 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* = 105 , a = 10-10 , so ist zwar ||r|| = ax* 10-5, trotzdem ist der Fehler || x - x* || = 105. 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.
Im obigen Beispiel erhält man für kleines a für k(A) 2 / a, d.h. das System ist für a << 1 extrem schlecht konditioniert. Anschaulich bedeutet das, daß die beiden Geraden g1: y = 0 und g2: y = ax nahezu parallel sind, sodaß auf Grund des ``schleifenden'' Schnitts von g1 und g2 eine kleine Änderung in den Ausgangsdaten eine große Änderung im Schnittpunkt (in der Lösung) nach sich zieht (Abb. 2).


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(Hn) wächst exponentiell mit n (Abb. 3). Zum Beispiel hat das System Hn ·x = b für bi = nj=1 hij 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.

hmat1cond.png
Abb. 3:  Konditionszahlen der Hilbertmatrizen Hn.




File translated from TEX by TTH, version 3.06.
On 24 May 2008, 21:50.