Lineare Gleichungssysteme mit n Unbekannten lösen

[d3g]

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]

[CalganX]

mschnell hat die Unit von d3g noch ein wenig verändert. Die entsprechende Unit befindet sich im Anhang, inkl. eines Demoprojektes.

Erläuterungen zur Unit:

Die Gauss/Matrix Unit bietet einige Funktionen zur Mathematik mit Matrizen dynamischer Größe. U.a. Lösung von linearen Gleichungssystemen nach Gauss-Algorithmus. Es werden auch mehrere Gleichungssysteme in einem Schritt berechnet.

Delphi-Quellcode:

// Copyright: Julian und Michel Schnell, Krefeld, Germany, [email][e-mail][/email]

interface
type
TVector = array of Extended;
TMatrix = array of TVector;

function SolveLinearSystem(A: TMatrix): TVector;
// sloves Ax=y, with A a square martrix
// Parameters:
// A: Matrix with the vector y added "at the right side" to the original A
// so A has a dimension of n, n+1
// Result: solution vector x

function SolveLinearSystems(A: TMatrix): TMatrix;
// sloves AX=Y, with A a square martrix, X and Y matrices with the same hight as A
// Parameters:
// A: Matrix with the matrix y added "at the right side" to the original A
// so A has a dimension of n, n+m (n = height of A, m = height of Y)
// Result: solution matrix X

function MatrixTrans (A: TMatrix) : TMatrix;
// transposes a matrix

function MatrixMulti (A, B: TMatrix) : TMatrix;
// multiplies matrices
// Parameters:
// A and B matrices, width of A = height of A
// Result: A*B

function MatrixVektorMulti (A: TMatrix; V:TVector ) : TVector;
// multiplies a matrix and a vector
// Parameters:
// A matrix, width of A = height of v
// Result: A*v

function MatrixConcat (A, B: TMatrix): TMatrix;
// Adds B to the right side of A
// Parameters:
// matrices A abd B with the same height

function MatrixIdent (n: Integer): TMatrix;
// returns an identity matrix I with dimension n

function MatrixInvert (A: TMatrix) : TMatrix;
// inverts a matrix
// Parameters:
// A a square matrix
// Result: A^-1

procedure MatrixBubbleSort(var A: TMatrix; m: Integer);
// sorts the lines of a matrix so that the column m is ascending
// (in place sort algorithm)