Programmumsetzungen

Die drei obigen Algorithmen per Hand auszuführen wäre sehr zeitaufwendig und mühselig. Ständig diesselben Rechenoperationen zu betreiben, wären für einen Menschen zu ermüdend. Ein einziger Rechenfehler würde die ganze Arbeit wertlos machen, da dieser Rechenfehler das Ergebnis völlig verfälschen würde.
Mehr an Bedeutung gewinnen die Algorithmen, wenn man sie von Computern ausführen läßt. So lassen sich Rechenfehler ausschließen, man muß lediglich noch mit Rundungsfehlern zurechtkommen. Außerdem ist die Ausführungszeit fast vernachlässigbar, in Bruchteilen von Sekunden wird das Ergebnis vom Computer berechnet.
Weil der Computer einen Algorithmus nicht einfach bearbeiten kann, muß man diesen erst in eine für den Computer verständliche Form umsetzten. Die Umsetzung ist dann ein Computerprogramm.
Die Entwicklungen für Computerprogramme der drei beschriebenen Algorithmen soll nachfolgend vollzogen werden. Die Programme werden in PASCAL geschrieben und sollten systemunabhängig sein.
Im folgenden Verlauf wird zunächst das Grundprinzip und dann erst das fertige Programm, das, wie das zugehörige Ablaufprotokoll, im Anhang zu finden ist. Es empfiehlt bei der Programmbeschreibung, den Programmquelltext vor sich liegen zu haben, um die einzelnen Schritte besser verfolgen zu können.
Bei der Programmierung wurde darauf geachtet, daß sich die einzelnen Programme sehr ähneln. Es wurde versucht ähnliche Programmteile wie z.B. Ein- und Ausgabe einheitlich zu gestalten, folglich werden gleiche Programmteile auch nur einmal im Text erklärt.


1. Gauß\'sches Eliminationsverfahren


1.1 Grundprinzip

Dem Computer muß zu Beginn erst einmal das Gleichungssystem übergeben werden, auf das das Gauß\'sche Eliminationsverfahren angwendet werden soll. Die Befehlsanweisungen sollte man dazu sinnvollerweise in einer Prozedur mit dem Namen \"Eingabe\" zusammenfassen. Jetzt ist es wieder nützlich, sich das Gleichungssystem in Matrix-Schreibweise vorzustellen. Lediglich Koeffizientenmatrix und den Freigliedervektor muß dem Computer übergeben werden. Es bietet sich an, diese Werte in einem zwei-dimensionalen nx(n+1) Array vom Typ Real abzulegen. In der Spalte n+1 legt man die Freiglieder ab. Diese Feld muß ebenso global definiert sein wie die Dimension n des Gleichungssystems.
Die Eingaberoutine könnte dann wie folgt lauten:



PROCEDURE Eingabe;
BEGIN
Write(\'Anzahl der Gleichungen: \');
ReadLn(n);
FOR i:=1 TO n DO
FOR j:=1 TO n+1 DO
BEGIN
Write(\'Koeffizient[\',i,\',\',j,\'] \');
ReadLn(Koeffizient[i,j]);
END;
END;
Dem Computer ist das Gleichungssystem nun bekannt, der Eliminationsvorgang kann gestartet werden. Die Prozedur \"Eliminieren\" soll die Befehlsanweisungen dafür in sich vereinen.
Insgesamt benötigt man drei Laufschleifen. Die äußerste Schleife, die k-Schleife, führt die k-Eliminationsstufen und Pivotierungen durch. Die k-Schleife läuft von eins bis n-1, da in einem n-dimensionalen Gleichungssystem n-1 Eliminationsstufen notwendig sind, um dieses in ein gestaffeltes Gleichungsystem zu überführen. Vor jeder Eliminationsstufe muß zunächst der Pivotierungsvorgang durchgeführt werden. Dieser ist unbedingt notwendig, da das Programm sonst eventuell mit \"Division by Zero\" abbrechen könnte.
Die Pivotierung gestaltet sich ziemlich einfach. Für das Pivotelement nimmt man zu Anfang Null an, sucht dann mit einer Schleife die k-te Spalte nach vom Betrag her größeren Werten ab und ermittelt dann die Zeile des größten Wertes.
Diese Zeile vertauscht man dann in einer anderen Laufschleife mit der k-ten Zeile. Ist das Pivotelement trotzdem null, so muß der Eliminationsvorgang abgebrochen werden, da das Gleichungssystem in diesem Fall keine eindeutige Lösung besitzt.
Erst jetzt kann eliminiert werden. Dazu benötigt man zwei verschachtelte Laufschleifen. Die äußere sei die i-Schleife, sie gibt an, in welcher Zeile gerade eliminiert wird. In dieser wird der Faktor c ermittelt, mit dem dann die Zeile k multipliziert und von der Zeile i subtrahiert wird. Dieser Vorgang wird in der innersten Schleife, der j-Schleife, ausgeführt.
Die Prozedur \"Eliminieren\" ist damit komplett:



PROCEDURE Eliminieren;
VAR Pivotelement, Faktor: REAL;
Pivotzeile, i, j, k: INTEGER;
BEGIN
FOR k:= 1 TO n - 1 DO {* Gesamtelimination}
BEGIN
Pivotelement:= 0;
Pivotzeile:= k;
FOR i:= k TO n DO {* Pivotelement suchen}
IF Abs(Koeffizient[i, k]) Pivotelement THEN
BEGIN
Pivotelement:= Abs(Koeffizient[i, k]);
Pivotzeile:= i;
END;
IF k Pivotzeile THEN {* wenn k=i, Pivotelement bereits richtig}
FOR j:= k TO n + 1 DO {* Zeilen vertauschen}
Exchange(Koeffizient[k, j], Koeffizient[Pivotzeile, j]);

IF Koeffizient[k, k] = 0 THEN Abbruch;
FOR i:= k + 1 TO n DO {* Elimination auf Stufe k}
BEGIN
Faktor:= Koeffizient[i, k] / Koeffizient[k, k];
FOR j:= k TO n + 1 DO
Koeffizient[i,j]:= Koeffizient[i,j]-Faktor*Koeffizient[k,j];
END;
END;
END;
Die Vorwärtselimination ist beendet, das Rücksubstituieren kann gestartet werden. Die dazugehörige Prozedur soll \"Rückwärtseinsetzen\" heißen. Man benötigt zwei Laufschleifen, eine äußere i- und eine innere j-Schleife. Die i-Schleife durchläuft rückwärts die Werte von n bis eins.
Um den folgenden Schritt besser verstehen zu können, ruft man sich am besten noch einmal das gestaffelte Gleichungssystem aus I.2.1 ins Gedächtnis:

a11x1 + a12x2 +...+ a1nxn = a1,n+1
a22x2 +...+ a2nxn = a2,n+1

..................................
annxn = an,n+1
Löst man nun jeweils die i-te Gleichung nach xi auf, so erhält man folgendes:
x1 = (a1,n+1 - (a12x2 +...+ a1nxn)):a11

x2 = (a2,n+1 - (a23x3 +...+ a2nxn)):a22

........................................

xn = (an,n+1 - 0):ann
Die j-Schleife berechnet nun für die Zeile i die Summe in der innersten Klammer (oben kursiv). Die j-Schleife durchläuft die Werte i+1 bis n.
Man könnte denken, es gäbe ein Problem, wenn i=n ist, weil dann j=n+1 sein muß und somit auch xn+1 herangezogen würde, das natürlich nicht existiert. Doch löst sich dieses Problem in Wohlgefallen auf, da die j-Schleife in diesem Fall überhaupt nicht mehr durchlaufen wird, da dann der Startwert j=n+1 größer als der Endwert n ist.
Im weiteren Verlauf der i-Schleife wird nun noch geprüft, ob ann=0 und gegebenenfalls wieder ein \"Division by Zero\" Fehler abgefangen.
Letztendlich wird die Lösung xi berechnet.
Als Programmtext sieht dies wie folgt aus:



PROCEDURE RückwärtsEinsetzen;
VAR i, j, NullPos: INTEGER;
Term: REAL;
BEGIN
FOR i:= n DOWNTO 1 DO
BEGIN
Term:= 0;
FOR j:= i + 1 TO n DO
Term:= Term + Koeffizient[i, j] * Lösung[j];

IF Koeffizient[i, i] = 0 THEN Abbruch;
Lösung[i] := (Koeffizient[i, n + 1] - Term) / Koeffizient[i, i];
END;
END;
Jetzt steht die simultane Lösung in dem Array \"Lösung\". Dieses kann man sich am Ende noch mit einer Prozedur \"DruckeLösung\" ausgeben lassen. Das komplette Programm liegt im Anhang als Quelltext vor.

1.2 Programmbeschreibung

Die Programmbeschreibung wird in die verschiedenen Unterpunkte Variablen und Konstanten, Prozeduren und Funktionen und Hauptprogramm unterteilt.
Das unten beschriebene Programm gibt zusätzlich auf Wunsch noch Rechenzwischenschritte aus, in der Programmbeschreibung wird darauf allerdings nicht eingegangen. Anzumerken ist noch, daß diese Zwischenschrittausgabe sehr einfach gehalten wurde, d.h. es kann zu Ausgaben kommen, die mathematisch nicht korrekt sind (z.b. zwei aufeinanderfolgende Verknüpfungszeichen).
Die einzelnen Schritte werden in der Reihenfolge beschrieben, in der sie auch im Programmtext zu finden sind.


Globale Variablen und Konstanten
In der Konstante MaxGleichungen wird die maximale Anzahl der Gleichungen, d.h. die maximale Dimension des Gleichungssystems festgelegt. Dies ist notwendig, da in PASCAL die Größe von Arrays nur über Konstanten definiert werden kann.
Die Konstante ZeichenBreite enthält die Breite des Bildschirms in Zeichen. Dieser Wert wird verwendet, um Daten zentriert ausgeben zu können.
Koeffizient ist eine Variable vom Typ Matrix in diesem zweidimensionalen Feld wird die Koeffizientenmatrix zusammen mit den Freigliedern gespeichert.
Lösung soll am Ende des Programms die simultane Lösung des Gleichungssystems enthalten. Lösung kann genauso wie Koeffizient Zahlen vom Typ Real speichern.
n beeinhaltet die tatsächliche Dimension des Gleichungssystems, allerdings darf n MaxGleichungen nicht überschreiten.

Prozeduren
Die Prozedur Eingabe liest, wie der Name bereits ausdrückt, die vom Benutzer gewünschten Daten ein und speichert diese dem entsprechend in n und Koeffizient.
DruckeGleichung druckt die Gleichung Nr auf den Bildschirm. Diese Prozedur wird von DruckeGleichungssystem verwendet, um das komplette Gleichungssystem in einem ansehnlichen Format auszugeben.
Die Prozedur Abbruch übernimmt die Fehlerausgabe, falls das eingegebene Gleichungssystem keine eindeutige Lösung besitzt.

Eliminieren und RückwärtsEinsetzen sind die Prozeduren, die vorher bereits ausführlich entwickelt wurden. Sie wurden jedoch noch um einige Zeilen erweitert, damit die jeweiligen Zwischenschritte mit ausgegeben werden können, sonst wurde nichts verändert.
DruckeLösung übernimmt die Ausgabe der simultanen Lösung.

Hauptprogramm:
Das Hauptprogramm konnte sehr kurz gehalten werden, da die überwiegende Arbeit bereits von den Prozeduren übernommen wird.
Zu Beginn des Hauptprogramms wird mittels der Prozedur Eingabe die Eingabe des Benutzers bearbeitet, darauf wird der Bildschirm gelöscht und die eingegebenen Koeffizienten als Gleichungssystem mit DruckeGleichungssystem dargestellt. Schließlich übernehmen die Prozeduren Eliminieren und RückwärtsEinsetzen das Lösen des Gleichungssystems. Diese Lösung wird dann zum Schluß mit DruckeLösung ausgegeben.

2. Gauß-Jordan Verfahren

2.1 Grundprinzip

Das PASCAL-Programm für das Gauß-Jordan-Verfahren läßt sich sehr gut aus dem gerade besprochenen entwickeln. Da sich der Gauß-Jordan Algorithmus vom Gauß\'schen Eliminationsvefahren erst ab dem Jordan\'schen Reduktionsvefahren unterscheidet, lassen sich die Prozeduren \"Eingabe\" und \"Eliminieren\" unverändert übernehmen.
Man muß lediglich eine \"GaussJordanReduktion\"-Prozedur programmieren. Dies gestaltet sich auch relativ einfach, da die \"Eliminieren\"-Prozedur dafür nur etwas umgeschrieben werden muß:



PROCEDURE GaussJordanReduktion;
VAR i, j, k: INTEGER;
Faktor: REAL;
BEGIN
FOR k:= n DOWNTO 1 DO
BEGIN
IF Koeffizient[k, k] = 0 THEN Abbruch;
FOR i:= k - 1 DOWNTO 1 DO
BEGIN
Faktor:= Koeffizient[i, k] / Koeffizient[k, k];
FOR j:= n + 1 DOWNTO i DO
Koeffizient[i,j]:=Koeffizient[i,j]-Faktor*Koeffizient[k, j];
END;
Lösung[k] := Koeffizient[k, n + 1] / Koeffizient[k, k];
END;
END;
Die Pivotierung kann völlig entfallen. In dem gestaffelten Gleichungssystem existiert pro Spalte jeweils nur ein Pivotelement, so daß es also gar keine Alternative zum Tauschen gibt. Totales Pivotieren wäre möglich, würde das Programm aber nur unnötig verlängern.
Die weiteren Änderungen sind ebenfalls sehr gering. Es müssen lediglich die Grenzen und die Laufrichtung der k-, i- und j-Schleife geändert werden (kursiv hervorgehoben).

2.2 Programmbeschreibung

Auch dieses Programm liegt wieder im Anhang als Source vor, wobei auch dieses Programm wieder Zwischenschritte ausgeben kann. Auf die Dokumentation dieser Programmteile wird wie oben verzichtet.
Die Programmbeschreibung kann sehr kurz gehalten werden, da dieser Quelltext mit dem Programm Gauß in den Punkten Variablen und Konstanten, Funktionen und Hauptprogramm fast übereinstimmt. Es wurde lediglich die Prozedur Rückwärtseinsetzen durch die Prozedur GaussJordanReduktion ausgetauscht, sonst blieb alles beim Alten.

3. Gauß-Seidel Verfahren

3.1 Grundprinzip

Das Programm für dieses Verfahren unterscheidet sich, genauso wie das Verfahren selbst, von den beiden vorherigen fast völlig. Folglich läßt sich von den beiden vorherigen, abgesehen von der Ein- und Ausgaberoutine und Pivotierung, nichts übernehmen.
Bevor das Iterationsverfahren starten kann, muß zuerst das beste Pivotelement gesucht werden, um ein divergieren des Iterationsprozesses zu vermeiden. Im Gegensatz zu den beiden vorherigen Verfahren wird beim Gauß-Seidel Verfahren in einem Zug durchpivotiert, d.h. die Pivotierung steht getrennt vor den Iterationsanweisungen.
Zur Pivotierung benötigt man drei Schleifen, wobei zwei in einer verschachtelt sind. Die äußerste Schleife ist die k- Schleife, die von 1 bis n-1 läuft und die Spalten durchläuft, aus welchen das Pivotelement ermittelt werden soll. Der nun folgende Programmteil aus i- und j-Schleife ist aus den beiden vorher beschriebenen Programmen bekannt.
Nun kommt die Überprüfung des Gleichungssystems. Die äußere i-Schleife durchläuft die zeilenweise das Gleichungssystem. Stößt die innere j-Schleife nur auf eine Zeile, in der das Pivotelement nicht überwiegt, ist das Gleichungssystem bereits ungeeignet, bricht die i-Schleife ab und eine Meldung wird ausgegeben. An dieser Stelle könnte man auch gleich das Programm abbrechen.
Die Variable SZeile enthält immer die Summe der Beträge der Koeffizienten der Zeile i ohne das Diagonalelement.
Aber nun zur Umsetzung des Iterationsprozesses. Zuerst einmal muß man die initiale Näherung bestimmen, d.h. alle Felder des Arrays \"Lösung\" auf Null setzen, dann kann man den Iterationsprozeß starten. Man benötigt zur Programmierung drei ineinander verschachtelte Schleifen. Die äußerste bestimmt die Dauer des Iterationsprozesses. Dieser soll nämlich, nach folgenden Kriterien abgebrochen werden. Entweder wenn die Lösung gefunden wurde, d.h. wenn sich die Näherung gegenüber der vorherigen nicht mehr verbessert, oder wenn die Zahl der Iterationen über einen vorgegebenen Wert hinaus gewachsen ist. Hierzu bietet sich eine REPEAT-DO-UNTIL-Schleife an, da man die Anzahl der Durchläufe vorher nicht kennt. Zu Beginn dieser Schleife wird die Laufvariable um
eins erhöht, dann startet die nächste Schleife, die i-Schleife. Diese beginnt bei dem Wert eins und endet bei n. In dieser Schleife werden jeweils die xi mit i=1..n berechnet. Jetzt folgt bereits die innerste Schleife, die j-Schleife, die ebenfalls die Werte von eins bis n durchläuft. Mit dieser Schleife wird jeweils der Zähler aus den vorherigen Iterationen bzw. der initialen Näherung und den dazugehörigen Koeffizienten berechnet. Nach dieser Schleife kann nun das neue xi berechnet werden, wobei vorher wieder ein \"Division by Zero\"-Fehler abgefangen werden muß. Nun wird noch verglichen, ob sich Genauigkeit der neuberechneten Lösung gegenüber der alten verbessert hat. Ist dies nicht der Fall, so wird der Iterationsvorgang beendet. Die andere Möglichkeit wäre, daß die Maximalzahl der Iterationen überschritten wurde und somit keine genaue Lösung gefunden wurde. Diesen Fall mit einzubeziehen ist wichtig, da das Programm vom Computer bis in die Ewigkeit weitergeführt würde, ohne eine Lösung zu finden. Hier der Programmtext:



PROCEDURE GaussSeidel;
VAR Pivotelement, Zähler, NeueLösung, SZeile: REAL;
Pivotzeile, i, j, k, Iteration: INTEGER;
Fertig, geeignet: BOOLEAN;
BEGIN
FOR k:= 1 TO n - 1 DO
BEGIN
Pivotelement:= 0;
FOR i:= 1 TO n DO {* Pivotelement suchen}
IF Abs(Koeffizient[i, k]) Pivotelement THEN BEGIN
Pivotelement:= Abs(Koeffizient[i, k]);
Pivotzeile:= i;
END;
IF k Pivotzeile THEN {* wenn k=i, Pivotelement bereits günstig}
FOR j:= 1 TO n + 1 DO {* Zeilen vertauschen}
Exchange(Koeffizient[k, j], Koeffizient[Pivotzeile, j]);
END;

geeignet:= TRUE; {* überprüfen, ob Gleichungssystem geeignet}
i:= 1;
REPEAT
SZeile:= 0.0;
FOR j:= 1 TO n DO
IF i j THEN SZeile:= SZeile + Abs(Koeffizient[i, j]);
IF SZeile = Abs(Koeffizient[i, i]) THEN geeignet:= FALSE;
i:= i + 1;
UNTIL (i n) OR NOT geeignet;
IF NOT geeignet THEN
WriteLn(\'Gleichungssystem ungeeignet\');
WriteLn;

FOR i:= 1 TO n DO
Lösung[i] := 0; {* erste Lösungsnäherung}

Iteration:= 0;

REPEAT
Iteration:= Iteration + 1;
Fertig:= TRUE;
FOR i:= 1 TO n DO
BEGIN
Zähler:= Koeffizient[i, n + 1];
FOR j:= 1 TO n DO
IF i j THEN Zähler:= Zähler - Koeffizient[i, j]*Lösung[j];
IF Koeffizient[i, i] = 0 THEN Abbruch;
NeueLösung:= Zähler / Koeffizient[i, i];
IF NeueLösung Lösung[i] THEN Fertig:= FALSE;
Lösung[i] := NeueLösung;
END;
UNTIL Fertig OR (Iteration = MaxIterationen);
END;

3.2 Programmbeschreibung

Auch das letzte Programm befindet sich wieder im Anhang. Beim Betrachten des Sourcecodes erkennt man wieder den bekannte Programmaufbau und die Struktur von Gauss. Allerdings ersetzt in diesem Fall die Prozedur GaussSeidel gleich die zwei Prozeduren Eliminieren und RückwärtsEinsetzen.
Außerdem gibt es noch eine Konstante MaxIterationen. Diese Konstante speichert die maximal zulässige Anzahl der Iterationen.
Der Rest ist mit den beiden vorherigen Programmtexten identisch.
3.3 Schlußbemerkung

Damit ist nun auch der letzte Teil der Facharbeit abgehandelt. Es wurde versucht die Algorithmen so knapp und genau wie möglich vorzustellen. Mit diesen Algorithmen sollten alle linearen NxN-Gleichungssysteme zu lösen sein.
Zuletzt soll noch geklärt werden, wann welches Verfahren eingesetzt werden soll.
Die Verfahren von Gauß und Jordan sind universell einsetzbar, d.h. sie finden immer eine Lösung. Besonders bei großen Gleichungssystemen steigt die Bearbeitungszeit allerdings rapide an, da beim Gauß\'schen Eliminationsverfahren die Anzahl Arbeitsschritte für ein Gleichungssytem der Dimension N mit N3/3 ansteigt. Das Gauß-Jordan Verfahren benötigt sogar noch mehr Arbeitschritte, es arbeitet ungefähr anderthalb mal langsamer als das Verfahren von Gauß. Das Gauß\'sche Verfahren ist eines der schnellsten Lösungsverfahren, es kann höchstens von iterativen, z.B. Gauß-Seidel (N2-3), übertroffen werden.
Das Gauß-Seidel-Verfahren ist bei sehr großen Gleichungssystemen meistens das schnellste, wenn es funktioniert.
Neben diesen Algorithmen existieren natürlich noch viele weitere, wie z.B. Cramer\'sche Regel oder LU-Dekomposition, die ihrerseits auch wieder Vor- und Nachteile haben.