Robert Sedgewick: Algorithmen

[ Inhaltsverzeichnis ] [ vorhergehende Seite ] [ nächste Seite ] [ Stichwort ]


37. Gaußsches Eliminationsverfahren



Beschreibung des Verfahrens

Im allgemeinen Fall möchten wir ein System von N Gleichungen mit N Unbekannten lösen:

a11x1 + a12x2 + ... + a1NxN = b1,
a21x1 + a22x2 + ... + a2NxN = b2,
.
.
aN1x1 + aN2x2 + ... + aNNxN = bN.

Diese Gleichungen lassen sich bei Benutzung der Matrixschreibweise als eine einzige Matrixgleichung schreiben:

Formel

oder einfach Ax = b, wobei A die Matrix darstellt, x die Variablen und b die rechten Seiten der Gleichungen. Da die Zeilen von A zusammen mit den Elementen von b umgeformt werden, ist es zweckmäßig, b als die (N+1)-te Spalte von A anzusehen und ein Feld vom Typ N mal (N+1) zu verwenden, um beides zu speichern.

Nunmehr kann die Etappe der Vorwärts-Elimination wie folgt zusammengefaßt werden: Zuerst eliminiere man durch Addition geeigneter Vielfacher der ersten Gleichung zu jeder der anderen Gleichungen die erste Variable in allen Gleichungen mit Ausnahme der ersten. Dann eliminiere man durch Addition geeigneter Vielfacher der zweiten Gleichung zu jeder der Gleichungen von der dritten bis zur N-ten die zweite Variable in allen Gleichungen mit Ausnahme der ersten zwei. Dann eliminiere man die dritte Variable in allen Gleichungen mit Ausnahme der ersten drei usw. Um die i-te Variable in der j-ten Gleichung (für j zwischen i + 1 und N) zu eliminieren, multiplizieren wir die i-te Gleichung mit aji/aii und subtrahieren sie von der j-ten Gleichung. Dieser Prozeß wird in kürzerer Form durch den folgenden Programmabschnitt beschrieben:

    for i:=1 to N do
      for j:=i+1 to N do
        for k:=N+1 downto i do
          a[j,k]:=a[j,k]-a[i,k]*a[j,i]/a[i,i];

Das Programm besteht aus drei ineinander verschachtelten Schleifen, und die Gesamtlaufzeit ist im wesentlichen proportional zu N3. Die dritte Schleife wird rückwärts durchlaufen, so daß vermieden wird, daß a[j,i] zerstört wird, bevor es benutzt wird, um die Werte anderer Elemente in der gleichen Zeile zu korrigieren.

Der Programmabschnitt im obigen Absatz ist zu einfach, um völlig richtig zu sein: a[i,i] könnte den Wert null haben, so daß eine Division durch null auftreten könnte. Dies läßt sich jedoch leicht vermeiden, da wir jede beliebige Zeile (von der (i+1)-ten bis zur N-ten) mit der i-ten Zeile vertauschen können, damit a[i,i] in der äußeren Schleife von null verschieden wird. Falls keine solche Zeile gefunden werden kann, ist die Matrix singulär: Sie besitzt keine eindeutige Lösung. (Unser Programm könnte dies explizit melden, oder wir könnten den Fehler eventuell als eine Division durch null anzeigen lassen.) Zu dem obigen Programmabschnitt müssen wir weitere Programmzeilen hinzufügen, um eine weiter unten befindliche Zeile zu finden, die in der i-ten Spalte ein von null verschiedenes Element aufweist, und um diese Zeile dann mit der i-ten Zeile zu vertauschen. Das Element a[i,i], das schließlich benutzt wird, um die von null verschiedenen Elemente unterhalb der Diagonale in der i-ten Spalte zu eliminieren, wird Pivotelement genannt.

In Wirklichkeit ist es ratsam, noch etwas mehr zu tun, als nur eine Zeile mit einem von null verschiedenen Wert in der i-ten Spalte zu finden. Es ist am besten, diejenige Zeile (von der (i + 1)-ten bis zur N-ten) zu verwenden, für die der Wert in der i-ten Spalte dem absoluten Betrag nach am größten ist. Der Grund dafür ist, daß bei der Berechnung erhebliche Fehler entstehen können, wenn der Pivotwert, der verwendet wird, um eine Zeile mit einem Faktor durchzumultiplizieren, sehr klein ist. Wenn a[i,i] sehr klein ist, wird der Faktor a[j,i]/a[i,i], der benutzt wird, um die i-te Variable aus der j-ten Gleichung (für j von i + 1 bis N) zu eliminieren, sehr groß. Tatsächlich kann der Fall eintreten, daß er so groß wird, daß die eigentlichen Koeffizienten a[j,k] ihm gegenüber derart klein werden, daß der Wert von a[j,k] durch »Rundungsfehler« verzerrt wird.

Einfach ausgedrückt, Zahlen, die sich hinsichtlich ihrer Größe beträchtlich unterscheiden, können in dem gewöhnlich zur Darstellung reeller Zahlen verwendeten System von Gleitkommazahlen nicht exakt addiert oder subtrahiert werden, und die Verwendung eines kleinen Pivotwertes erhöht die Wahrscheinlichkeit beträchtlich, daß solche Operationen ausgeführt werden müssen. Die Verwendung des größten Wertes in der i-ten Spalte von der (i+1)-ten bis N-ten Zeile gewährleistet, daß der Faktor, mit dem die Zeile durchmultipliziert wird, stets kleiner als 1 ist, und beugt dieser Art von Fehlern vor. Man könnte erwägen, außerhalb der i-ten Spalte zu suchen, um ein großes Element zu finden, doch es ist bewiesen worden, daß exakte Lösungen erhalten werden können, ohne daß man diese zusätzliche Komplikation in Kauf nehmen muß.

Das folgende Programm für die Etappe der Vorwärts-Elimination des Gaußschen Eliminationsverfahrens ist eine einfache Implementation dieses Prozesses. Für jedes i von 1 bis N durchlaufen wir die i-te Spalte, um das größte Element zu finden (in den Zeilen nach der i-ten Zeile). Die Zeile, die dieses Element enthält, wird mit der i-ten Zeile vertauscht, und danach wird die i-te Variable in den Gleichungen i+1 bis N genau wie zuvor eliminiert:

    procedure eliminate;
      var i,j,k,max: integer;
          t: real;
      begin
      for i:=1 to N do
        begin
        max:=i;
        for j:=i+1 to N do
          if abs(a[j,i])>abs(a[max,i]) then max:=j;
        for k:=i to N+1 do
          begin t:=a[i,k]; a[i,k]:=a[max,k]; a[max,k]:=t end;
        for j:=i+1 to N do
          for k:=N+1 downto i do
            a[j,k]:=a[j,k]-a[i,k]*a[j,i]/a[i,i];
        end
      end;

In manchen Algorithmen wird gefordert, daß das Pivotelement a[i,i] benutzt wird, um die i-te Variable aus jeder Gleichung mit Ausnahme der i-ten (nicht nur aus der (i + 1)-ten bis N-ten) zu eliminieren. Dieser Prozeß wird vollständiges Pivotieren genannt; bei der Vorwärts-Elimination führen wir nur einen Teil dieser Operationen aus, daher wird der Prozeß partielles Pivotieren genannt. In Kapitel 43 betrachten wir einen Algorithmus für das vollständige Pivotieren.

Nachdem die Etappe der Vorwärts-Elimination abgeschlossen ist, enthält das Feld a unterhalb der Diagonale nur Nullen, und die Etappe des Rückwärts-Einsetzens kann ausgeführt werden. Das Programm hierfür ist sogar noch einfacher:

    procedure substitute;
      var j,k: integer;
      t: real;
      begin
      for j:=N downto 1 do
        begin
        t:=0.0;
        for k:=j+1 to N do t:=t+a[j,k]*x[k];
        x[j]:=(a[j,N+1]-t)/a[j,j]
        end
      end;

Ein Aufruf von eliminate, gefolgt von einem Aufruf von substitute, bewirkt die Berechnung der Lösung in dem aus N Elementen bestehenden Feld x. Bei singulären Matrizen könnte immer noch eine Division durch null auftreten; eine Bibliotheksroutine würde dies explizit prüfen. In Wirklichkeit führen die meisten Unterprogramme in Bibliotheken wesentlich umfangreichere Überprüfungen durch, worauf wir weiter unten noch eingehen werden.

Eigenschaft 37.1 Ein Gleichungssystem mit N Gleichungen und N Unbekannten kann unter Verwendung von ungefähr N3/3 Multiplikationen und Additionen gelöst werden.

Die Laufzeit von substitute ist O(N2), so daß die meiste Arbeit in eliminate ausgeführt wird. Eine Betrachtung dieser Routine zeigt, daß für jeden Wert von i die k-Schleife N - i + 2 mal und die j-Schleife N - i mal durchlaufen wird; das bedeutet, daß die innere Schleife Summe1<=i<=N(N - i + 2)(N - i) = N3/3 + O(N2) mal durchlaufen wird. Der Wert von -a[j,i]/a[i,i] kann außerhalb der k-Schleife berechnet werden, so daß die innere Schleife aus einer Multiplikation und einer Addition besteht. w.z.b.w.

Nachdem die Vorwärts-Elimination unter der Diagonale überall Nullen erzeugt hat, besteht eine andere mögliche Vorgehensweise darin, genau die gleiche Methode anzuwenden, um oberhalb der Diagonale überall Nullen zu erzeugen: Zuerst erzeuge man in der letzten Spalte überall Nullen, mit Ausnahme von a[N,N], indem man geeignete Vielfache von a[N,N] addiert, dann führe man das gleiche für die vorletzte Spalte aus usw. Das heißt, wir nehmen nochmals ein »partielles Pivotieren« vor, doch für den anderen »Teil« jeder Spalte, wobei wir die Folge der Spalten rückwärts durchlaufen. Nachdem dieser Prozeß (der Gauß-Jordansches Reduktionsverfahren genannt wird) abgeschlossen ist, sind nur noch Elemente auf der Diagonale von null verschieden, womit die Lösung unmittelbar gegeben ist. Allerdings ist bei diesem Prozeß die Anzahl der ausgeführten Rechenoperationen wesentlich größer als beim Rückwärts-Einsetzen.

Die Möglichkeit von Fehlern bei der numerischen Berechnung erfordert beim Gaußschen Eliminationsverfahren besondere Beachtung. Wie bereits erwähnt wurde, ist in Situationen Vorsicht geboten, in denen sich die Größen der Koeffizienten erheblich unterscheiden. Die Verwendung des größten zur Verfügung stehenden Elements in der Spalte für das partielle Pivotieren gewährleistet, daß beim Prozeß des Pivotierens nicht willkürlich große Koeffizienten erzeugt werden, doch es ist nicht immer möglich, größere Fehler zu vermeiden. Zum Beispiel treten sehr kleine Koeffizienten auf, wenn zwei verschiedene Gleichungen Koeffizienten haben, die sehr nahe beieinander liegen. In Wirklichkeit ist es allerdings möglich, im voraus zu bestimmen, ob solche Probleme zu ungenauen Werten in der Lösung führen werden. Jeder Matrix entspricht ein mit ihr verknüpfter Zahlenwert, der die Konditionszahl genannt wird und der benutzt werden kann, um die Genauigkeit der berechneten Lösung zu beurteilen. Eine gute Bibliotheksroutine für das Gaußsche Eliminationsverfahren berechnet neben der Lösung auch die Konditionszahl der Matrix, so daß die Genauigkeit der Lösung bestimmt werden kann. Eine vollständige Behandlung der damit zusammenhängenden Fragen würde über den Rahmen dieses Buches hinausgehen.

Beim Gaußschen Eliminationsverfahren mit partiellem Pivotieren unter Verwendung des größten zur Verfügung stehenden Pivotwertes kann »garantiert« werden, daß die erhaltenen Ergebnisse mit sehr kleinen Rechenfehlern behaftet sind. Es liegen sehr sorgfältig ausgearbeitete mathematische Ergebnisse vor, die zeigen, daß die berechnete Lösung sehr genau ist, außer für schlecht konditionierte Matrizen (was eher ein Hinweis auf Probleme im Gleichungssystem sein dürfte als im Lösungsverfahren). Der Algorithmus war Gegenstand sehr gründlicher theoretischer Untersuchungen und kann als ein numerisches Verfahren mit sehr breiten Anwendungsmöglichkeiten empfohlen werden.


[ Inhaltsverzeichnis ] [ vorhergehende Seite ] [ nächste Seite ] [ Stichwort ]