Lösungsverfahren zur numerischen Lösung linearer Gleichungssysteme

TGM_98 : LINGL.*, GAUSS.SBR, GAUSS.SYS

Ronald Hasenberger / cand.ing. TU-Wien

1. Allgemeines

Bei dem hier vorliegenden Text behandle ich zunächst möglichst kurz (und dadurch auch sehr ungenau) die Theorie der Lösung linearer Gleichungssysteme, um mich anschließend auf konkrete Lösungsverfahren, insbesondere das Gaußsche Eliminationsverfahren zu konzentrieren.

1.1. Notation

In diesem Text werden Matrizen durch Unterstreichen des jeweiligen Symbols gekennzeichnet z.B.:

A = Matrix A

Die Unterscheidung zwischen eindimensionalen Matrizen (Vektoren) und zweidimensionalen Matrizen (Matrizen im eigentlichen Sinn) erfolgt dadurch, daß Vektoren durch Kleinbuchstaben bezeichnet werden.

Z.B.:

v = Vektor v

Soll auf einzelne Elemente der Matrix hingewiesen werden, erfolgt das in Indexschreibweise mit Kleinbuchstaben z.B.:

┌ ┐

│ a11 a12 a13 ... a1n │

│ a21 a22 a23 ... a2n │

A = │ a31 a32 a33 ... a3n │

│ : : : : │

│ am1 am2 am3 ... amn │

└ ┘

Die hier vorgestellte Matrix ist eine m x n Matrix mit m Zeilen und n Spalten.

Eine derartige Matrix kann man sich nun wieder aufgespalten in ihre Zeilen und Spalten vorstellen, wobei hierdurch m Zeilen- und n Spaltenvektoren definiert werden. Als Spaltenvektor wird im folgenden ein Vektor bezeichnet, der nur aus einer einzigen Spalte (m x 1 Matrix) besteht z.B.:

┌ ┐

│ s1 │

│ s2 │

s = │ s3 │

│ : │

│ sm │

└ ┘

Als Zeilenvektor wird im folgenden ein Vektor bezeichnet, der nur aus einer einzigen Zeile (1 x n Matrix) besteht z.B.:

┌ ┐

z = │ z1 z2 z3 ... zn │

└ ┘

Ein lineares Gleichungssystem läßt sich formal als die Matrizenmultiplikation einer Koeffizientenmatrix A mit einem Spaltenvektor x, welche als Ergebnis wieder einen Spaltenvektor b ergibt, darstellen:

A * x = b

ausführlich:

a11 x1 + a12 x2 + a13 x3 + ... + a1n xn = b1

a21 x1 + a22 x2 + a23 x3 + ... + a2n xn = b2

a31 x1 + a32 x2 + a33 x3 + ... + a3n xn = b3

: : : : :

am1 x1 + am2 x2 + am3 x3 + ... + amn xn = bm

Das hierdurch erklärte System ist durchaus auch für nicht quadratische Koeffizientenmatrizen erklärt, wobei dann aber für die Lösbarkeit bzw. eindeutige Lösbarkeit bestimmte Voraussetzungen getroffen werden müssen, auf die ich hier nicht näher eingehen möchte (Näheres siehe (1)). Im folgenden werde ich mich ausschließlich auf quadratische Koeffizientenmatrizen (d.h. m = n) beziehen.

1.2 Lösbarkeit von linearen Gleichungssystemen

Damit ein derartiges Gleichungssystem eindeutig lösbar ist (was uns letztlich für diesen Bereich hauptsächlich interessiert), müssen einige Voraussetzungen betreffend die Koeffizientenmatrix festgehalten werden.

1.2.1 Lineare Abhängigkeit von Vektoren

Sind a1, a2, a3, ..., an beliebige Vektoren, so wird die mit beliebigen Zahlen µ1, µ2, µ3, µn gebildete Summe

µ1 a1 + µ2 a2 + µ3 a3 + ... + µn an

eine Linearkombination der Vektoren ak genannt. Diese Linearkombination wird im n-dimensionalen Raum durch einen Polygonzug repräsentiert.

Man kann sich nun fragen, ob es durch geeignete Wahl der Konstanten µk möglich ist, einen geschlossenen Polygonzug zu erzeugen d.h.:

µ1 a1 + µ2 a2 + µ3 a3 + ... + µn an = 0

Dies ist offensichtlich möglich, indem alle µk zu Null gesetzt werden, wobei dann aber auch der Polygonzug (zu einem Punkt) entartet.

Ist es aber auch für µk möglich, die nicht alle gleich Null sind, so spricht man von linearer Abhängigkeit zwischen den Vektoren ak. In diesem Fall läßt sich zumindest einer der Vektoren durch eine Linearkombination der restlichen Vektoren darstellen.

Z.B. sind die Vektoren

┌ ┐ ┌ ┐ ┌ ┐

│ 1 │ │ 4 │ │ 7 │

a1 = │ 2 │; a2 = │ 5 │; a3 = │ 8 │

│ 3 │ │ 6 │ │ 9 │

└ ┘ └ ┘ └ ┘

linear abhängig, weil die Linearkombination

µ1 a1 + µ2 a2 + µ3 a3 + ... + µn an

für µ1 = 1, µ2 = -2 und µ3 = 1 den Nullvektor ergibt, womit die Bedingung für lineare Abhängigkeit erfüllt ist.

Lineare Abhängigkeit ist immer dann gegeben, wenn eine Zeile (oder auch eine Spalte) der Koeffizientematrix ein Vielfaches einer anderen Zeile (Spalte) ist.

Die lineare Unabhängigkeit der Koeffizientenmatrix ist dann gegeben, wenn die Determinante der Koeffizientenmatrix ungleich Null ist.

1.2.2 Lösbarkeit von linearen Gleichungssystemen

Kehren wir nun wieder zurück zur Lösbarkeit von linearen Gleichungssystemen.

Es läßt sich hier zeigen, daß ein lineares Gleichungssystem nur dann eindeutig lösbar ist, wenn alle (n) Spalten- und Zeilenvektoren der Koeffizientenmatrix A linear unabhängig sind. Im Fall der Nichterfüllung dieser Bedingung ist das Gleichungssystem entweder nur nicht eindeutig lösbar (d.h. das Ergebnis ist wiederum ein geforderter Zusammenhang zwischen verschiedenen Variablen) oder nicht lösbar. Dies hängt hauptsächlich von der rechten Seite, dem Ergebnisvektor b ab.

2. Lösungsverfahren bei linearen Gleichungssystemen

2.1. Cramer'sche Regel

Die Cramer'sche Regel führt die Lösung eines linearen Gleichungssystemes auf die Berechnung von Determinanten zurück:

│ a11 ... a1 i-1 b1 a1 i+1 ... a1n │

1 │ a21 ... a2 i-1 b2 a2 i+1 ... a2n │

xi = ───── │ : : : : : │

Det A │ : : : : : │

│ an1 ... an i-1 bn an i+1 ... ann │

Hierbei wird zur Berechnung der i-ten Variable jeweils die i-te Spalte der Koeffizientenmatrix durch die linke Seite ersetzt, und die Determinante der so erhaltenen Matrix durch die Determinante der ursprünglichen Koeffizientenmatrix dividiert.

Aus dieser Form ist auch (numerisch) leicht ersichtlich, daß die Determinante der Koeffizientenmatrix nicht Null sein darf, da hier sonst eine Division durch Null erfolgen würde.

Aufgrund des sehr hohen Aufwands der Berechnung von Determinanten, und der Notwendigkeit der Berechnung von n + 1 Determinanten für die Lösung eines Gleichungssystems mit n Variablen wird dieses Verfahren selten angewendet, und wird hier auch nicht weiter behandelt.

2.2 Gaußsches Eliminationsverfahren

Das Gaußsche Eliminationsverfahren beruht darauf, daß ein Gleichungssystem (nebenbei: wie auch die Determinante einer Matrix) seinen Wert nicht ändert, wenn zu (von) einer Zeile ein Vielfaches einer anderen Zeile addiert (subtrahiert) wird.

Dies wird beim Gaußschen Eliminationsverfahren dazu benützt der Reihe nach die einzelnen Variablen zu eliminieren, indem man eine Gleichung als Hauptgleichung benützt, und zu jeder der restlichen Gleichungen ein geeignetes Vielfaches dieser Gleichung addiert, sodaß in diesen eine Variable nicht mehr vorkommt.

Wir betrachten hierzu das ursprüngliche System,

a11 x1 + a12 x2 + a13 x3 + ... + a1n xn = b1

a21 x1 + a22 x2 + a23 x3 + ... + a2n xn = b2

: : : : :

an1 x1 + an2 x2 + an3 x3 + ... + ann xn = bn

welches sich nach dem ersten Schritt der Transformation mit der ersten Gleichung als Hauptgleichung zu

a11 x1 + a12 x2 + a13 x3 + ... + a1n xn = b1

a22'x2 + a23'x3 + ... + a2n'xn = b2'

: : : :

an2'xn + an3'xn + ... + ann'xn = bn'

ergibt.

Diese Transformation ist in dieser Form nur für a11 <> 0 möglich, was aber keine ernsthafte Beschränkung darstellt, da durch Vertauschen zweier Zeilen, wodurch der Wert des Gleichungssystems ebenfalls nicht geändert wird, sicher eine Situation hergestellt werden kann, in der a11' <> 0 gilt; die Faktoren ergeben sich hierbei zu:

- a21 / a11 für die 2. Gleichung

- a31 / a11 für die 3. Gleichung

- an1 / a11 für die n-te Gleichung

Das selbe Verfahren kann man nun (für a22' <> 0) mit der 2. Gleichung als Hauptgleichung zur Elimination der Variablen x2 aus den Gleichungen 3 bis n verwenden.

Durch weitere Anwendung dieses Verfahrens kann man auf diese Art, wenn das System überhaupt eindeutig lösbar ist, das Gleichungssystem in die Form

r11 x1 + r12 x2 + r13 x3 + ... + r1n xn = c1

r22 x2 + r23 x3 + ... + r2n xn = c2

: : :

rnn xn = cn

überführen, aus dem man durch Lösen der letzten Gleichung xn (= cn / rnn), mit dem Einsetzen von xn in die n-1 - te Gleichung xn-1 und so weiter bis x1 erhalten kann.

Ist die Überführung des Gleichungssystems in die obengenannte Form nicht möglich, so ist das Gleichungssystem nicht eindeutig lösbar, d.h. singulär. Dieser Fall soll aus den weiteren Betrachtungen ausgeklammert werden, da dessen Behandlung einen wesentlich höheren Aufwand erfordern würde.

Man erhält dadurch ein neues, dem ursprünglichen Gleichungssystem aber äquivalentes Gleichungssystem, für welches ich folgende Notation benützen möchte:

R * x = c

Dieses Verfahren ist gleichwertig mit der Multiplikation der Koeffizientenmatrix A sowie der rechten Seite b mit einer unteren Dreiecksmatrix L von links (beachte: die Matrizenmultiplikation ist nicht kommutativ!). Die untere Dreiecksmatrix (eine Matrix, bei der sämtliche Elemente oberhalb der Hauptdiagonale gleich Null sind) ergibt sich dabei aus der Multiplikation mehrerer (n) unterer Dreiecksmatrizen:

L = Ln * ... * L2 * L1

Die Matrizen Lm haben dabei in der Hauptdiagonale nur Einsen stehen. Die Elemente unterhalb der Hauptdiagonale sind alle gleich Null, außer jenen in der m-ten Spalte, wo die Koeffizienten der Elimination

- am+1m / amm

- am+2m / amm

- an m / amm

stehen.

Dadurch ergibt sich z.B. die Matrix L1 zu:

┌ ┐

│ 1 0 0 ... 0 │

│ -a21 / a11 1 0 ... 0 │

│ -a31 / a11 0 1 ... 0 │

L1 = │ -a41 / a11 0 0 ... 0 │

│ : : : : │

│ -an1 / a11 0 0 ... 1 │

└ ┘

Die Matrizen R und c ergeben sich hierbei durch:

R = L * A

c = L * b

Diese Darstellung ist deshalb sehr günstig, weil sie die Berechnung eines Gleichungssystems für eine neue rechte Seite ermöglicht, ohne daß dafür die gesamte Elimination nochmals durchgeführt werden muß. Ich werde später darauf noch zurückkommen.

3. Numerik des Gaußschen Eliminationsverfahrens

Das Gaußsche Eliminationsverfahren, dessen Grundlagen im bisherigen Text bereits behandelt wurden ist grundsätzlich sehr einfach zu programmieren. Man findet dabei für eine einfache Version durchaus mit der Anwendung des bis jetzt Gesagten sein Auslangen.

Für höhere Ansprüche kann man das Verfahren aber noch verbessern.

Zunächst besteht hier die Möglichkeit die Genauigkeit des Verfahrens zu erhöhen:

3.1 Pivotisierung des Gleichungssystems

Hierzu möchte ich zunächst kurz auf die Zahlendarstellung im Computer eingehen.

Bei modernen Programmiersprachen wird in der Regel zwischen Integer (= ganzzahligen) und Real (= gebrochenen) Zahlen unterschieden.

Die Speicherung von Integerzahlen führt in der Regel kaum zu Problemen, da man sich hierbei auf einen Zahlenbereich beschränkt, der durch die verwendete Wortlänge exakt beschrieben werden kann.

Z.B. Integer = 16 Bit Worte ==>

Zahlenbereich: -2^(n-1) = -32768 <= Z <= 32767 = 2^(n-1) - 1)

Hierbei wird jede mögliche Zahl des Zahlenbereiches durch genau eine Bitkombination beschrieben, wodurch keine Fehler auftreten können.

Anders ist die Situation bei Realzahlen: Hier muß ein Dezimalbruch dargestellt werden, der im Extremfall nie endet (z.B. ã, 1/3, ...). Derartige Zahlen können nicht exakt dargestellt werden, sodaß man zur Verarbeitung im Computer nur die ersten n Stellen des Dezimalbruches verwendet (n ist gleich der Anzahl der gültigen Stellen). Um trotzdem große und kleine Zahlen darstellen zu können verwendet man hier eine Darstellung mit Mantisse und Exponent. Die Mantisse sind die eigentlichen Stellen der Zahl, der Stellenwert der gesamten Zahl ergibt sich durch Multiplikation mit 10^Exponent.

z.B.: Mantisse = 0,1234; Exponent = 2: Zahl = 0,1234 * 10^2 = 12,34

Eine Addition wird dabei so durchgeführt, daß die beiden Zahlen zunächst durch Verschieben der Stellen der kleineren Zahl auf den Exponenten der größeren Zahl gebracht und anschließend addiert werden. Hierbei kann es aber vorkommen, daß die eine Zahl soweit verschoben werden muß, daß keine gültige Stelle mehr für die Addition verfügbar ist, sodaß das Ergebnis so aussieht, als wenn Null addiert worden wäre.

z.B.:

a) Addition von 0,1234 * 10^0 + 0,5678 * 10^2 in einem Computer mit 5 gültigen Stellen; die führende Null wird hierbei nicht mitgezählt:

1. Schritt: Auf gleichen Exponenten bringen:

Hierbei werden beide Zahlen auf den Exponenten der größeren Zahl gebracht.

==> 0,12340 * 10^0 -> 0,00123 * 10^2

0,56780 * 10^2 -> 0,5678 * 10^2

Hierbei ging bereits die letzte gültige Ziffer des ersten Summanden verloren.

2. Schritt: Ziffernweise Addition

==> 0,00123 * 10^2

+ 0,56780 * 10^2

──────────────

0,56903 * 10^2

b) Addition von 0,1234 * 10^0 + 0,5678 * 10^8 in einem Computer mit 5 gültigen Stellen; die führende Null wird hierbei nicht mitgezählt:

1. Schritt: Auf gleichen Exponenten bringen:

Hierbei werden beide Zahlen auf den Exponenten der größeren Zahl gebracht.

==> 0,1234 * 10^0 -> 0,00000 * 10^8

0,5678 * 10^8 -> 0,56780 * 10^8

Hierbei ging die gesamte Information des ersten Summanden veloren.

2. Schritt: Ziffernweise Addition

Dieser Schritt ist in diesem Fall eine Addition von 0 zu 0,5678*10^8

Wie bereits in der Theorie erläutert, wird beim Gaußschen Eliminationsverfahren das Aij/Aii - fache einer Zeile (mit i: Nummer der zu eliminierenden Variablen, j: Nummer der Zeile, in der eliminiert werden soll) von der j - ten Zeile abgezogen. Wird dieser Faktor sehr groß, so "überdeckt" die Hauptgleichung die anderen Gleichungen. Im diesem Fall wird die Gleichung, in der eliminiert werden soll für den Rechner mit seiner begrenzten Darstellungsmöglichkeit nur mehr ein Vielfaches der Hauptgleichung sein, was die eindeutige Lösung des Gleichungssystems verhindern würde.

Dieser Faktor wird nun (betragsmäßig) genau dann sehr groß, wenn Aii betragsmäßig sehr klein wird.

Dies ist auch die Möglichkeit zu einer Abhilfe: Wenn als Aii das (betragsmäßig) größte aller Aij mit j von i bis Anzahl der Gleichungen gewählt wird, werden die Faktoren klein gehalten (betragsmäßig kleiner oder gleich 1).

Die Elimination der Variablen wird daher so ablaufen, daß man vor der Elimination einer Variablen das betragsmäßig größte Element der Koeffizientenmatrix in der Spalte der betrachteten Variablen unter und auf der Hauptdiagonale sucht, diese Gleichung zur Hauptgleichung erklärt, und mit Hilfe dieser Gleichung die Variable in den anderen Gleichungen eliminiert.

Für die Verarbeitung im Rechner, bei der man die zu eliminierenden Variablen in der Regel in Form einer Schleife von 1 bis Anzahl der Variablen eliminieren wird, ergibt sich daher die Notwendigkeit, das Gleichungssystem umzustellen.

Hier besteht die überlegungsmäßig einfachste Möglichkeit die i-te Zeile mit der Zeile, in der das betragsmäßig größte Element gefunden worden ist, auszutauschen, indem die entsprechenden Elemente der Koeffizientenmatrix A und des Spaltenvektors b kopiert werden. Aufgrund der hohen Anzahl von Zuordnungsvorgängen ist diese Methode aber relativ langsam.

Eine zweite Möglichkeit ist die Verwendung eines Steuervektors.

3.1.1 Pivotisierung mittels Steuervektor

Hierbei bedient man sich folgender Struktur:

s1 -> a11 a12 .... a1n b1

s2 -> a21 a22 .... a2n b2

: : : : :

sn -> an1 an2 .... ann bn

Man definiert einen Steuervektor s, der ebensoviele Zeilen hat wie die Koeffizientenmatrix A. Bei der Adressierung der Zeilen der Koeffizientenmatrix A sowie des Spaltenvektors b adressiert man nun die zeilen nicht direkt, sondern über den Steuervektor:

z.B. A(.s(.i.),j.), wenn man die das j-te Element der i-ten Zeile adressieren will; ebenso beim Spaltenvektor b:

b(.s(.i.).), wenn man die i-te Zeile des Spaltenvektors b adressieren will.

Vor Beginn der Elimination ist dafür zu sorgen, daß das i-te Element des Steuervektors den Wert i enthält, wodurch gewährleistet ist, daß wirklich die erste Zeile adressiert wird, wenn A(.s(.1.),j.) adressiert wird.

Wenn nun aber im Zuge einer Pivotsuche die gewünschte Hauptgleichung z.B. die 5. Gleichung statt der 3. ist, so reicht ein Vertauschen von s(.5.) mit s(.3.) um den Zeilentausch durchzuführen. Es darf trotzdem nicht direkt, wie in diesem Beispiel naheliegend zu sein scheint

s(.5.) := 3 und

s(.3.) := 5

durchgeführt werden, da es ja möglich ist, daß durch eine frühere Vertauschung s(.5.) gar nicht mehr 5 enthält (desgleichen mit s(.3.))!

3.2 Verbesserung des Ergebnisses

Jener Ergebnisvektor x, den man nach der Durchführung des Gaußschen Eliminationsverfahrens und anschließender Berechnung der einzelnen Koordinaten erhält ist zwangsläufig nicht genau.

Es besteht jetzt die Möglichkeit, durch Einsetzen des Ergebnisvektors x in die ursprüngliche (!!!!) Gleichung einen Fehler für die linke Seite zu bestimmen.

diff = b - A * xe

xe: errechneter Ergebnisvektor

Diese Differenz ist nun aber noch kein Maß für den Fehler der Ergebnisvektors. Aufgrund der Distributivitätseigenschaften der Matrizenmultiplikation gilt aber:

A * x - A * xe = b - A * xe = diff

woraus folgt:

A * (x - xe) = diff oder

A * Delta_x = diff

x: exakter Ergebnisvektor

Delta_x: Fehler, um den der errechnete Ergebnisvektor zu klein ist.

Daraus folgt nun aber, daß durch Lösung des Gleichungssystems

A * Delta_x = diff

eine (wiederum nicht exakte) Lösung für den Fehler berechnet werden kann, der bei der Lösung des Gleichungssystems begangen wurde.

Indem dieser Fehler zur ursprünglichen Lösung addiert wird, erhält man eine genauere Lösung des ursprünglichen Systems. Diese Iterationen sind aber nur in einem beschränkten Ausmaß sinnvoll, da ja die Genauigkeit auch durch die Faktoren, mit denen die Elimination durchgeführt wurde beeinflußt wird, deren Genauigkeit aber durch diese Methode nicht gesteigert werden kann.

Zur Berechnung des Fehlervektors Delta_x gibt es nun die einfache, aber aufwendige, Möglichkeit die gesamte Elimination nochmals durchzuführen.

Schneller ist es, wenn man unter Ausnützung einiger Eigenschaften der Elimination die Berechnung für eine neue linke Seite einfacher durchführt. Der Schlüssel hierzu liegt in der Möglichkeit, die Elimination als Multiplikation sowohl der Koeffizientenmatrix A als auch des Spaltenvektors b mit einer unteren Dreiecksmatrix L von links zu interpretieren:

R = L * A = Ln * ... * L2 * L1 * A

c = L * b = Ln * ... * L2 * L1 * b

wodurch sich das Gleichungssystem, wie bei der Theorie bereits erläutert als

R * x = c

darstellt.

Die Matrix R haben wir hierbei bereits einmal berechnet, es ist dies eine obere Dreiecksmatrix, die sich bei der Elimination von A ergibt.

Dadurch reicht es, die Matrizenmultiplikation von L mit b und anschließend die Berechnung der xi aus dem reduziertzen System durchzuführen.

Dieses Matrizenmultiplikation kann mit der Kenntnis aller Lm ebenso gut wie mit der Kenntnis der unteren Dreiecksmatrix L durchgeführt werden. Wenn mann alle Matrizen Lm in einem Feld speichert ergibt sich eine untere Dreiecksmatrix:

┌ ┐

│ 1 0 0 ... 0 │

│ -a21 / a11 1 0 ... 0 │

│ -a31 / a11 -a32'/ a21' 1 ... 0 │

│ -a41 / a11 -a42'/ a21' -a43"/ a31" ... 0 │

│ : : : : │

│ -an1 / a11 -an2'/ a21' -an3"/ a31" ... 1 │

└ ┘

Die "interessanten" Bereiche (jene, in denen eine nicht von vornherein bekannte Information steht) dieser Matrix sind angenehmerweise genau jene, die bei der reduzierten Koeffizientenmatrix Null werden, sodaß jene Speicherplätze der Koeffizientenmatrix, die zu Null werden, zur Speicherung eben jener Teile dieser Matrix verwendet werden können.

In etwas weniger mathematischer Form kann man auch sagen, daß man die Quotienten, mit denen man eine Variable aus einer Gleichung eliminiert hat, am Platz des Koeffizienten für diese Variable in der entsprechenden Gleichung abspeichert, da der Koeffizient nach der Elimination Null wird, was aufgrund der Struktur des Lösungsverfahrens bekannt ist.

Hierdurch ergibt sich nach Durchführung der gesamten Reduktion mit Speicherung der Quotienten folgender Zustand der Matrix:

┌ ┐

│ a11 a12 a13 ... a1n │

│ a21 /a11 a22' a23' ... a2n' │

A'= │ a31 /a11 a32'/a2' a33" ... a3n" │

│ : : : ... : │

│ an1 /a11 an2'/a22' an3"/a33" ... ann"' │

└ ┘

Bei der Berechnung der cn wird nun sinnvollerweise wieder das Feld für die bn doppelt verwendet, wodurch sich als eine Möglichkeit zur Berechnung der cn folgende Methode ergibt, wobei nicht mehr zwischen bn und cn unterschieden wird, da nur ein Array vorhanden ist.

┌──────────────────────────────────────────────┐

│ I := 1 to Dim - 1 │

│ ┌────────────────────────────────────────────┤

│ │ J := I + 1 to Dim │

│ │ ┌──────────────────────────────────────────┤

│ │ │ b (.J.) := b (.J.) - A (.J,I.) * b (.I.) │

└─┴─┴──────────────────────────────────────────┘

4. Weitere Anwendungen des Gaußschen Eliminationsverfahrens

4.1 Berechnung von Determinanten

Für die Berechnung von Determinanten mit Hilfe von Computern ist der wesentlichste Satz jener, daß bei einer oberen oder unteren Dreiecksmatrix der Wert der Determinante gleich dem Produkt der Elemente der Hauptdiagonale ist.

Das Gaußsche Eliminationsverfahren ist nun ein Verfahren, welches eine Matrix in eine obere Dreiecksmatrix umwandelt, sodaß (im wesentlichen) nur mehr die Elemente der Hauptdiagonale miteinander multipliziert werden müssen.

Das Vorzeichen der Determinante ist nach dieser Maßnahme aber noch nicht unbedingt korrekt, da eine Determinante bei einer Vertauschung von Zeilen oder Spalten ihr Vorzeichen wechselt, das heißt, es ist notwendig das Produkt der Elemente der Hauptdiagonale mit (-1)^Anzahl der Vertauschungen zu multiplizieren.

4.2 Inversion von Matrizen

Matrizeninversion kann am einfachsten durchgeführt werden, indem man n Gleichungssysteme in n Variablen löst:

A * X = E

X ist hierbei eine n x n Matrix, die Inverse zu A. Die linke Seite ist in diesem Fall ebenfalls eine Matrix, die Einheitsmatrix, bei der alle Elemente außer den Elementen der Hauptdiagonale 0 sind; letztere sind 1.

┌ ┐

│ 1 0 0 ... 0 │

│ 0 1 0 ... 0 │

E = │ 0 0 1 ... 0 │

│ : : : : │

│ 0 0 0 ... 1 │

└ ┘

Diese Gleichungssysteme kann man auf 2 Arten relativ effizient lösen:

Die erste Möglichkeit wäre die spaltenweise Bestimmung der Inversen X indem in ein Gleichungssystem

A * x = b

der Reihe nach als rechte Seite die Spalten der Einheitsmatrix eingesetzt werden. Die m-te Spalte der Inversen ist hierbei die Lösung des obigen Gleichungssystems mit der m-ten Spalte der Einheitsmatrix (alle Elemente gleich 0; nur das m-te gleich 1) als rechter Seite. Hierbei müßte nur einmal das gesamte Eliminationsverfahren abgearbeitet werden, danach könnte man das Verfahren zur Berechnung eines Gleichungssystems mit geänderter rechter Seite verwenden.

Dieses Verfahren bietet sich insbesondere an, wenn man auf ein vorhandenes System von Unterprogrammen zur Durchführung des Gaußschen Eliminationsverfahrens zurückgreifen kann.

Die zweite Möglichkeit ist wahrscheinlich etwas schneller, benötigt aber doch einige spezielle Unterroutinen:

Hierbei wird das gesamte System

A * X = E

betrachtet. Multipliziert man dieses System von links mit einer unteren Dreiecksmatrix L vom selben Aufbau wie beim Gaußschen Eliminationsverfahren, so erhält man folgendes System

R * X = L,

wobei die Matrix R die in eine obere Dreiecksmatrix umgerechnete Koeffizientenmatrix ist. Hierbei kann man die Spalten der Inversen X erhalten, indem man die Rückrechnung (!) des Gaußschen Eliminationsverfahrens für ein Gleichungssystem durchführt, dessen rechte Seite die entsprechende Spalte der unteren Dreiecksmatrix L ist. Die Dreiecksmatrix L müßte hierbei aber explizit berechnet werden.

5. Hinweise zur Programmierung

Es empfiehlt sich bei der Programmierung eines derartigen Programmes, eine Routine zu schreiben, die einem die Gestalt des aktuellen Gleichungssystems zeigt.

Im Falle der Verwendung der Version mit Steuervektor bedeutet dies, daß sowohl die Information vorhanden sein sollte, welche Zeile des ursprünglichen Systems eine ausgedruckte Zeile ist, als auch welche Zeile des aktuellen Systems sie darstellt, oder, was die selbe Information beinhaltet, wo die Quotienten enden und die eigentliche Koeffizientenmatrix beginnt.

Bei der Durchführung des Eliminationsverfahrens ist weiter darauf zu achten, daß eine Division durch Null verhindert (abgefangen) wird, wobei die Abfrage besser nicht auf Null erfolgen sollte, sondern eine, vom Problem abhängige Schwelle gesetzt werden sollte, unter der angenommen wird, daß dieses Element in Wirklichkeit (bei unendlich genauer Rechnung) Null wäre. In diesem Fall ist das Gleichungssystem nicht eindeutig lösbar.

6. Quellennachweis

(1) H.J. Dirschmid; Mathematische Grundlagen der Elektrotechnik

2. Auflage 1987;

Verlag Vieweg

(2) Mathematik A2 für Elektrotechnik

Vorlesung 101.695 Sommersemester 1987/88 Prof. Dirschmid

Technische Universität Wien

(3) Programmiertechnik für Elektrotechniker

Vorlesung 356.895 Wintesemester 1988/89 Prof. Barth

Technische Universität Wien

GAUSS.SYS

(* Dieses File definiert alle nötigen Typen für das Programm Gauss *)

Const MaxN = 20; (* Maximalanzahl der Gleichungen und Variablen *)

Minimum = 1E-30; (* Grenze für die Erkennung einer singulären *)

(* Matrix *)

Type Koeff = Array (.1..MaxN, 1..MaxN.) of Real;

(* Typ für die Speicherung der Koeffizienten *)

Spalte = Array (.1..MaxN.) of Real;

(* Typ zur Speicherung der rechten Seite und *)

(* des Ergebnisvektors *)

ISpalte = Array (.1..MaxN.) of Integer;

(* Typ zur Speicherung des Steuervektors *)

GAUSS.SBR

(* Dieses Unterprogrammpaket berechnet ein lineares Gleichungssystem

der Form

Ax = B,

A ... Matrix n x n (n: Dimension des Systems)

B ... Ergebnisvektor (Matrix n x 1)

x ... Variablenvektor (Matrix n x 1) *)

Procedure Gauss_Druck (DA: Koeff; DB:Spalte; DSteuer: ISpalte; DDimension: Integer);

(* Procedure zum Drucken eines gestaffelten Systems, wobei zwischen den

Quotienten und der eigentlichen Matrix ein ▌ gezeichnet wird *)

Var I, J: Integer; (* Zählervariablen *)

Dummy: Char; (* Dummy Zeichen zum Warten auf die Fortsetzung *)

Begin (* Gauss_Druck *)

Writeln;

For I := 1 to DDimension do

Begin For J := 1 to DDimension do

Begin If DSteuer (.I.) = J then Write ('▌ ')

else Write (' ');

Write (DA (.I,J.):7:1);

End;

(* Schreiben der Quotienten und der gestaffelten Matrix *)

Writeln (' = ',DB (.I.):7:1);

(* Schreiben des entsprechenden Elements, rechte Seite *)

End;

Writeln;

Writeln ('Weiter mit Return');

Read (Dummy);

(* nach dem Einlesen Fortsetzung des Programms *)

End; (* Gauss_Druck *)

Procedure Gauss_Init (Var ISteuer: ISpalte; IDimension: Integer);

(* Initialisierung des Steuervektors *)

Var I: Integer; (* Zählervariable *)

Begin (* Gauss_init *)

For I := 1 to IDimension do

ISteuer (.I.) := I;

End; (* Gauss_Init *)

Procedure Gauss_neue_Rechte (LA:Koeff; Var LB: Spalte; LSteuer: ISpalte; LDimension: Integer);

(* Diese Procedure berechnet die zur gestaffelten Matrix passende

rechte Seite aus der ursprünglichen rechten Seite des Systems

Var I, J: Integer; (* Zählervariablen *)

Begin (* Gauss_neue_Rechte *)

For I := 1 to LDimension - 1 do

For J := I + 1 to LDimension do

LB (.LSteuer (.J.).) :=

LB (.LSteuer (.J.).) - LA (.LSteuer (.J.),I.)

* LB (.LSteuer (.I.).);

End; (* Gauss_neue_Rechte *)

Procedure Gauss_Elim (Var GA:Koeff; Var GB: Spalte; Var GSteuer: ISpalte; GDimension: Integer; Var Nicht_Singulaer: Boolean);

(* Diese Procedure rechnet die Koeffizientenmatrix GA um in die

gestaffelte Koeffizientematrix, wobei der Steuervektor verändert

wird, und in den zu 0 gemachten Koeffizienten der Quotient Q

gespeichert wird *)

Var I, J: Integer; (* Zählervariablen *)

Max_Betrag: Real; (* maximaler Betrag, *)

(* der von Abs_Max gefunden wurde *)

Function Abs_Max (GAA: Koeff; Var GASteuer: ISpalte; GADimension: Integer; GAZeile: Integer): Real;

(* Diese Function bestimmt das betragsmäßig größte Element

der Koeffizientenmatrix GAA in der Spalte GAZeile unter der

Hauptiagonale unter Beachtung des Steuervektors GASteuer,

und tauscht die Zeile, in der das betragsmäßig größte Element

gefunden wurde durch Veränderung des Steuervektors in die

Zeile GAZeile. Der Rückgabewert von Abs_Max ist

der größte gefundene Betrag *)

Var I: Integer; (* Zählervariable *)

Hilf: Integer; (* Hilfsvariable zum Zeilentausch *)

Max: Real; (* größtes bisher gefundenes Element *)

Max_Zeile: Integer; (* Zeile, in der Max gefunden wurde *)

Begin (* Abs_Max *)

Max := ABS (GAA (.GASteuer (.GAZeile.),GAZeile.));

Max_Zeile := GAZeile;

For I := GAZeile + 1 to GADimension do

If ABS (GAA (.GASteuer (.I.),GAZeile.)) > Max then

Begin Max := ABS (GAA (.GASteuer (.I.),GAZeile.));

Max_Zeile := I;

End; (* Suche Zeile mit absolut größtem Element *)

If Max_Zeile <> GAZeile then (* vertausche die Zeilen *)

(* GAZeile und Max_Zeile *)

Begin Hilf := GASteuer (.GAZeile.);

GASteuer (.GAZeile.) := GASteuer (.Max_Zeile.);

GASteuer (.Max_Zeile.) := Hilf;

End;

Abs_Max := Max;

End; (* Abs_Max *)

Procedure Gauss_Single (Var GSA: Koeff; GSSteuer: ISpalte; GSDimension: Integer; GSZeile: Integer);

(* Diese Procedure macht einen Schritt einer Gauß-Transformation, d.h. sie

transformiert die Koeffizientenamtrix derart, daß in der Spalte GSZeile

unter der Hauptdiagonale nur mehr Nullen stehen würden.

Tatsächlich werden in jene Elemente aber die bei der Berechnung ver-

wendeten Quotienten geschrieben

Var I, J: Integer; (* Zählervariablen *)

Q: Real; (* Quotient *)

Begin (* Gauss_Single *)

For I := GSZeile + 1 to GSDimension do

(* I zählt die Zeilen *)

Begin Q := GSA (.GSSteuer (.I.),GSZeile.) /

GSA (.GSSteuer (.GSZeile.),GSZeile.);

GSA (.GSSteuer (.I.),GSZeile.) := Q;

(* Koeffizientenmatrix wird gleichzeitig

zum Speichern der Quotienten verwendet *)

For J := GSZeile + 1 to GSDimension do

(* J zählt die Spalten bei der Umrechnung einer Zeile

in die gestaffelte Form *)

GSA (.GSSteuer (.I.),J.) :=

GSA (.GSSteuer (.I.),J.) -

Q * GSA (.GSSteuer (.GSZeile.),J.);

End;

End; (* Gauss_Single *)

Begin (* Gauss_Elim *)

Nicht_Singulaer := true;

For I := 1 to GDimension do

If Nicht_Singulaer then

Begin Max_Betrag := Abs_Max (GA,GSteuer, GDimension, I);

(* suche abolut größtes Element in der Spalte I

unterhalb der Hauptdiagonale *)

(* und tausche die dazugehörige Zeile in die Zeile I *)

If not (Max_Betrag < Minimum) then

Gauss_Single (GA,GSteuer, GDimension, I)

(* führe einen Schritt der Transformation aus

else Nicht_Singulaer := false;

End;

Gauss_neue_Rechte (GA, GB, GSteuer, GDimension);

End; (* Gauss_Elim *)

Procedure Gauss_Rueck (RA: Koeff; RB: Spalte; RSteuer: ISpalte; RDimension: Integer; Var OK: Boolean; Var RX: Spalte);

(* Diese Procedure erledigt die Rückrechnung von einer gestaffelten Matrix

mit passender linker Seite auf den Variablenvektor x *)

Var I, J: Integer;

Hilf: Real; (* Hilfsvariable zum Berechnen des Divisors *)

Begin (* Gauss_Rueck *)

OK := true;

For I := RDimension downto 1 do

(* I zählt die Zeilen = Elemente des Ergebnisvektors *)

Begin Hilf := RB (.RSteuer (.I.).);

For J := RDimension downto I + 1 do

(* J zählt die Spalten bei der

Rücktransformation eines Elements

des Ergebnisvektors *)

Hilf := Hilf - RX (.J.) * RA (.RSteuer (.I.),J.);

If (Abs (RA (.RSteuer (.I.),I.)) > Minimum) then

RX (.I.) := Hilf / RA (.RSteuer (.I.),I.)

(* Verhinderung einer Division durch 0 *)

else OK := false;

End;

End; (* Gauss_Rueck *)

Procedure Gauss_Defekt (RkA: Koeff; RkB: Spalte; RkSteuer: ISpalte; RkDimension: Integer; RkX: Spalte; Var RkD: Spalte);

(* Gauss_Defekt berechnet die Differenz zwischen dem Vektor B

der Angabe und dem Vektor B, der bei einem Einsetzen

des berechneten Variablenvektors in die ursprüngliche Matrix A

entstehen würde *)

Var I, J: Integer; (* Zählervariablen *)

Begin (* Gauss_Defekt *)

For I := 1 to RkDimension do

Begin RkD (.I.) := RkB (.I.);

For J := 1 to RkDimension do

RkD (.I.) := RkD (.I.) - RkA (.I,J.) * RkX (.J.);

End;

End; (* Gauss_Defekt *)

LINGL.PAS

Program Lingl; (* Lösen eines linearen Gleichungssystems *)

(*$I GAUSS.SYS *) (* Definition der Typen für Gauss *)

Var Koeffmat: Koeff; (* Koeffizientenmatrix für Gauss *)

Rechte: Spalte; (* rechte Seite für Gauss *)

Koeffmat2: Koeff; (* original Koeffizientenmatrix

für Defektberechnung *)

Rechte2: Spalte; (* original rechte Seite

für Defektberechnung *)

Steuer: ISpalte; (* Steuervektor *)

Dim: Integer; (* Dimension des Gleichungssystems *)

X: Spalte; (* Ergebnisvektor des Gleichungssystems *)

Differenz: Spalte; (* Defekt des Gleichungssystems *)

DeltaX: Spalte; (* Fehler des Ergebnisvektors *)

OK: Boolean; (* Hilfsvariable, die anzeigt,

ob die letzten Operationen

ordnungsgemäß abliefen *)

I: Integer; (* Zählervariable zum allgemeinen Gebrauch*)

Dummy: Char; (* Dummychar zum Bremsen des Programms *)

(*$I GAUSS.SBR *) (* Prozeduren und Funktionen für Gauss *)

Procedure Eingabe (Var Matrix: Koeff;Var LSpalte: Spalte; Var Nbr: Integer; Var EingabeOK: Boolean);

(* Procedure zum Einlesen des Gleichungssystemes

Die Eingabe kann wahlweise von einem Textfile oder der Tastatur erfolgen.

Das Textfile muß dabei alle Informationene beinhalten, die sonst vom Benutzer

abgefragt würden:

o) Dimension des Gleichungssystems

o) Koeffizientenmatrix

o) rechte Seite

Aufbau des Textfiles:

Dimension

a11 a12 b1

a21 a22 b2 *)

Var Eing_File: Text; (* File, von dem gelesen werden kann *)

Eing_Filenam: String (.11.);

(* Name des Eingabefiles *)

I, J: Integer; (* Zählervariablen *)

Eing_Quelle: Char; (* Single Char für div. Abfragen *)

E_File: Boolean; (* true, wenn vom File gelesen werden soll *)

Begin (* Eingabe *)

EingabeOK := true;

Write ('Eingabe von File oder von Tastatur (F/T): ');

Readln (Eing_Quelle);

If (Eing_Quelle in (.'F','f'.)) then E_File := true

else E_File := false;

If E_File then Begin Repeat

Write ('Eingabefilename: ');

Readln (Eing_Filenam);

Writeln;

Assign (Eing_File, Eing_Filenam);

(*$I-*) Reset (Eing_File);(*$I+*)

Until (IOResult = 0);

Readln (Eing_File,Dim);

If Dim > MaxN then

Begin

Writeln ('Dimension des Gleichungssystems

zu groß');

EingabeOK := false;

Delay (10000);

End;

End

else Begin Repeat Write

('Dimension des Gleichungssystems:');

Readln (Nbr);

Writeln;

Until Nbr <= MaxN;

End;

For I := 1 to Dim do

Begin For J := 1 to Dim do

If E_File then Read (Eing_File,Matrix (.I,J.))

else Begin Write ('A [',I:1,',',J:1,'] = ');

Readln (Matrix (.I,J.));

End;

If E_File then Readln (Eing_File, LSpalte (.I.))

else Begin Write ('B [',I:1,'] = ');

Readln (LSpalte (.I.));

Writeln;

End;

If E_File then Close (Eing_File);

End; (* Eingabe *)

Begin (* Lingl *)

ClrScr;

Eingabe (Koeffmat, Rechte, Dim, OK);

ClrScr;

If OK then

Begin Koeffmat2 := Koeffmat;

Rechte2 := Rechte;

(* sichern der ursprünglichen Daten für spätere *)

(* Verwendung *)

Gauss_Init (Steuer,Dim);

Writeln ('Eingegebenes Gleichungssystem');

Gauss_Druck (Koeffmat, Rechte, Steuer, Dim);

Writeln;

(* Ausgabe der ursprünglichen Matrix *)

Gauss_Elim (Koeffmat, Rechte, Steuer, Dim, OK);

(* berechnen eines ersten Ergebnisses *)

If OK then

(* nur, wenn Elimination ordnungsgemäß *)

Begin

Writeln ('Gleichungssystem nach der Elimination:');

Gauss_Druck (Koeffmat, Rechte, Steuer, Dim);

(* Rücktransformation des Systems *)

Gauss_Rueck

(Koeffmat, Rechte, Steuer, Dim, OK, X);

Writeln;

Writeln

('Ergebnis unmittelbar nach Gleichungslösung');

Writeln;

Writeln ('Ergebnisvektor:');

Writeln;

For I := 1 to Dim do

Writeln ('X',I:1,' = ',X (.I.));

(* Berechnung des Defekts mit ursprünglicher *)

(* Koeffizientenmatrix *)

Gauss_Defekt

(Koeffmat2, Rechte2, Steuer, Dim, X, Differenz);

(* Ausgabe des Defekts *)

Writeln;

Writeln ('Defekt:');

Writeln;

For I := 1 to Dim do

Writeln ('Defekt',I:1,' = ',Differenz (.I.));

(* Aufbau eines neuen Gleichungssystems mit Fehlern *)

(* als rechter Seite und Berechnung desselben *)

(* Elimination nicht mehr nötig, da sich die *)

(* Koeffizientenmatrix nicht geändert hat *)

Gauss_neue_Rechte (Koeffmat, Differenz, Steuer, Dim);

(* Ruecktransformation -> DeltaX *)

Gauss_Rueck

(Koeffmat, Differenz, Steuer, Dim, OK, DeltaX);

Writeln;

Writeln ('Fehler des Ergebnisvektors:');

Writeln;

For I := 1 to Dim do

Begin

Writeln ('DeltaX',I:1,' = ',DeltaX (.I.));

X(.I.) := X(.I.) + DeltaX (.I.);

End;

Writeln;

Writeln ('Weiter mit Return');

Read (Dummy);

Writeln ('Ergebnis nach 1 Iteration:');

Writeln;

For I := 1 to Dim do Writeln ('X',I:1,' = ',X (.I.));

Gauss_Defekt

(Koeffmat2, Rechte2, Steuer, Dim, X, Differenz);

Writeln;

Writeln ('Defekt:');

Writeln;

For I := 1 to Dim do

Writeln ('Defekt',I:1,' = ',Differenz (.I.));

Gauss_neue_Rechte (Koeffmat, Differenz, Steuer, Dim);

Gauss_Rueck

(Koeffmat, Differenz, Steuer, Dim, OK, DeltaX);

Writeln;

Writeln ('Fehler des Ergebnisvektors');

Writeln;

For I := 1 to Dim do

Writeln ('DeltaX',I:1,' = ',DeltaX (.I.));

Writeln;

Writeln ('Weiter mit Return');

Read (Dummy);

End

else Begin Writeln;

Writeln ('Fehler: singuläre Matrix');

Writeln;

Writeln ('Berechnung abgebrochen');

End;

End. (* Lingl *)