Lösen eines linearen Gleichungssystems mit n Unbekannten
1. Wozu braucht man das?
Keine Ahnung *g*. Ich hab's einfach als Herausforderung gesehen einen Algorithmus zu fidnen, der diese Aufgabe absolviert und den dann auf Object Pascal umzusetzen. Vielleicht macht einer von euch das ja auch stangenweise an der Schule/FH/Uni und kann es deshalb brauchen...
Zu dieser Beschreibung muss ich noch sagen, dass ich das noch nicht in der Schule gehabt habe und dass ich deshalb Fehler gemacht haben könnte, schreibt mir bitte, wenn ihr einen findet. Ich habe das ganze Verfahren aber natürlich durchgetestet und der Code funktioniert auch.
2. Definitionen
Da ich aus eigener Erfahrung weiß, dass es für Uneingelesene schwer ist einen mathematischen Text zu lesen, weil man andauernd nochmal nachschauen muss "War das schon definiert?" / "Was war das nochmal?". Deshalb hier alle Definitionen gleich zu Anfang.
Code:
---------------------------------------------------------------------------------
| A | eine m-mal-n-Matrix (m Zeilen, n Spalten), die aus |
| | dem linearenGleichungssystem entsteht |
---------------------------------------------------------------------------------
| A' | die schon durch die sogennante Vorwärtselimination |
| | abgeänderte Matrix |
---------------------------------------------------------------------------------
| m, n | Spalten- und Zeilenanzahl der Matrizen A und A' |
---------------------------------------------------------------------------------
| a(1; 1), a(1; 2), | die Zellen der Matrix A (bis auf die letzte Spalte) |
| ..., a(m; n - 1) | |
---------------------------------------------------------------------------------
| a'(1; 1), a'(1; 2), | die Zellen der schon abgeänderten Matrix A' (bis auf |
| ..., a'(m; n - 1) | die letzte Spalte) |
---------------------------------------------------------------------------------
| c(1), c(2), ..., c(m) | die Werte der einzelnen Gleichungen des Systems auf |
| | der anderen Seite des Gleichheitszeichens, diese |
| | Werte sind auch in der Matrix A in der letzten |
| | Spalte vorhanden, also ist c(a) = a(a, n) |
---------------------------------------------------------------------------------
| c'(1), c'(2), ..., | dasselbe wie oben nur für die Matrix A' |
| c'(m) | |
---------------------------------------------------------------------------------
| x(1), x(2), ..., x(m) | die Unbekannten des Gleichungsystems |
---------------------------------------------------------------------------------
| i, k | Zählvariablen |
---------------------------------------------------------------------------------
3. Der Gaußsche Algorithmus
Der Gaußsche Algorithmus ist der Algorithmus, der hier verwendet wird, um die Gleichunssysteme zu lösen. Dazu wird das Gleichungssystem in eine Matrix übertragen. Als Beispiel für diese Bescheibung werde ich das folgende Gleichunssystem verwenden:
Code:
I. x(1) + 2*x(2) - x(3) = -8
II. 2*x(1) - x(2) + x(3) = 11
III. 2*x(1) - 2*x(2) - 2*x(3) = 2
Also ersten Schritt trägt man in diese Matrix A die Koeffizienten der Variablen des Gleichungssystems sowie die Lösungen ein, sodass die Matrix nachher so aussieht:
Code:
{ a(1; 1) a(1; 2) ... a(1; n - 1) c(1)
A = . . . .
. . . .
a(m; 1) a(m; 2) ... a(m; n - 1) c(m) }
Für unser Beispiel ist A also:
Code:
{ 1 2 -1 -8
A = 2 -1 1 11
2 -2 -2 2 }
3.1 Vorwärtseliminination
Um das Gleichungssystem nach dem Gaußschen Algorithmus zu lösen, muss man aus A eine linke obere Matrix A' machen (eine Matrix, bei der alle Elemente unterhalb der Hauptdiagonalen 0 sind), dabei darf das Gleichungssystem natürlich nicht so abgeändert werden, dass nachher x(1) bis x(m) andere Werte haben.
Dies wird im ersten Schritt des Gaußschen Algorithmus getan:
1. Eliminationsschritt:
Man subtrahiere das a(2; 1) / a(1; 1)-fache der ersten Zeile von der zweiten Zeile, das a(3; 1) / a(1; 1)-fache der ersten Zeile von der dritten Zeile usw.
Unsere Matrix sieht nach diesem ersten Eliminationsschritt wie folgt aus:
Code:
{ 1 2 -1 -8
A' = 0 -5 3 27
0 -6 0 18 }
Wie man sieht, entehen in der ersten Spalte (außer bei a'(1; 1)) Nullen, damit kommen wir unserem Ziel einer linken oberen Matrix A' näher.
2. Eliminationsschritt:
Das Verfahren ist das gleiche wie beim ersten Schritt, jedoch benutzt man als Multiplikator für die 2., 3., ..., m. Zeile nun nicht a(2; 1) / a(1; 1) etc. sondern a'(3; 2) / a'(2; 2), a'(4; 2) / a'(2; 2) etc. Nach diesem Schritt sieht die Matrix nun wie folgt aus:
Code:
{ 1 2 -1 -8
A' = 0 -5 3 27
0 0 -3,6 -14,4 }
Für unsere Matrix ist kein weiterer Eliminationsschritt nälog, für größere Matrizen wird einfach nach dem obigen System bis zum m - 1. Eliminatioinsschritt fortgefahren.
3.2 Pivotisierung
Das obere Verfahren scheitert in einem Fall: wenn ein Element der Hauptdiagonalen a'(i; i) = 0 ist (Division by Zero). Deshalb muss man vor jedem Eliminationsschritt prüfen, ob nicht dieser Fall eintritt. Wenn ja, dann sucht man sich unter den noch folgenden Zeilen einfach die Zeile k heraus, bei der |a'(k; i)| am größten ist. Ist kein Element a(k; i) <> 0 vorhanden, ist das Gleichungssystem unlösbar.
3.3 Rückwärtssubstitution
Wenn man die unter Punkt 3.1 gewonnene Matrix A' wieder in ein Gleichungssystem umwandelt, erkennt man sehr schnell, wozu dieser Schritt gut war:
Code:
I. x(1) + 2*x(2) - x(3) = -8
II. -5*x(2) + 3*x(3) = 27
III. -3,6*x(3) = -14,4
Man sieht: das Gleichungssystem ist so gut wie gelöst, man muss nur noch durch die Koeffizienten dividieren und die vorherigen x einsetzen.
Das ist schnell in eine Funktion umgesetzt:
Code:
n
____
\
x(i) = (c'(i) / a'(i; i)) - > x(k) * (a'(i; k) / a'(i; i))
/___
k = i + 1
Daraus ergibt sich:
x(1) = 2
x(2) = -3
x(3) = 4
4. Umsetzung nach Object Pascal
Dies Nach Object Pascal umzusetzen ist keine große Schwierigkeit mehr.
Der folgende Code zeigt die universelle Funktion und wie sie für unser Beispiel aufgerufen wird:
Delphi-Quellcode:
program gauss;
{$APPTYPE CONSOLE}
uses
SysUtils;
type
TGaussSolved = array of Extended;
TGaussLine = TGaussSolved;
TGaussMatrix = array of TGaussLine;
function SolveLinearSystem(A: TGaussMatrix; m, n: Integer): TGaussSolved;
var
i, j, k: Integer;
Pivot: TGaussLine;
PivotRow: Integer;
Multiplicator, Sum: Extended;
begin
SetLength(A, m, n);
for i := 0 to m - 1 do
// Vorwärtselimination
for j := i to m - 2 do begin
if (A[j, j] = 0) then begin
// Pivotisierung
SetLength(Pivot, n + 1);
Pivot := A[j];
PivotRow := 0;
for k := j + 1 to m - 1 do begin
if (Abs(A[k, j]) > Abs(Pivot[j])) then begin
Pivot := A[k];
PivotRow := k;
end;
if (PivotRow > 0) then begin
A[PivotRow] := A[j];
A[j] := Pivot;
end else
raise EMathError.Create('System insolvable');
end;
end;
Multiplicator := A[j + 1, i] / A[i, i];
for k := i to n - 1 do
A[j + 1, k] := A[j + 1, k] - (Multiplicator * A[i, k]);
end;
// Rückwärtssubstitution
SetLength(Result, m);
for i := m - 1 downto 0 do begin
Sum := 0;
Result[i] := 0;
for k := i to m - 1 do
Sum := Sum + Result[k] * A[i, k] / A[i, i];
Result[i] := A[i, n - 1] / A[i, i] - Sum;
end;
end;
var
A: TGaussMatrix;
Res: TGaussSolved;
i, j: Integer;
begin
SetLength(A, 3, 4);
A[0][0] := 1; A[0][1] := 2; A[0][2] := -1; A[0][3] := -8;
A[1][0] := 2; A[1][1] := -1; A[1][2] := 1; A[1][3] := 11;
A[2][0] := 2; A[2][1] := -2; A[2][2] := -2; A[2][3] := 2;
for i := 0 to High(A) do begin
for j := 0 to High(A[i]) - 1 do
Write(FloatToStr(A[i, j]), '*x(', j + 1, ') + ');
WriteLn(#8#8, '= c(', i + 1, ')');
end;
Res := SolveLinearSystem(A, 3, 4);
WriteLn;
for i := 0 to High(Res) do
WriteLn('x(', i + 1, ') = ', FloatToStr(Res[i]));
ReadLn;
end.
Dann noch viel Spaß beim Lösen
,
d3g
[edit=sakura] Korrektur im Code. Mfg, sakura[/edit]