AGB  ·  Datenschutz  ·  Impressum  







Anmelden
Nützliche Links
Registrieren
Zurück Delphi-PRAXiS Code-Bibliothek Library: Algorithmen Delphi Lineare Gleichungssysteme mit n Unbekannten lösen

Lineare Gleichungssysteme mit n Unbekannten lösen

Ein Thema von d3g · begonnen am 22. Jun 2002 · letzter Beitrag vom 10. Jun 2007
Antwort Antwort
Benutzerbild von d3g
d3g

Registriert seit: 21. Jun 2002
602 Beiträge
 
#1

Lineare Gleichungssysteme mit n Unbekannten lösen

  Alt 22. Jun 2002, 09:18
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]
-- Crucifixion?
-- Yes.
-- Good. Out of the door, line on the left, one cross each.
  Mit Zitat antworten Zitat
CalganX

Registriert seit: 21. Jul 2002
Ort: Bonn
5.403 Beiträge
 
Turbo Delphi für Win32
 
#2

Re: Lineare Gleichungssysteme mit n Unbekannten lösen

  Alt 10. Jun 2007, 10:18
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]mschnell@bschnell.de[/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)
Angehängte Dateien
Dateityp: zip matrix_gauss_981.zip (6,9 KB, 397x aufgerufen)
»Mein Kaffee ist so schwarz — der fängt gleich an zu rappen...«
  Mit Zitat antworten Zitat
Themen-Optionen Thema durchsuchen
Thema durchsuchen:

Erweiterte Suche
Ansicht

Forumregeln

Es ist dir nicht erlaubt, neue Themen zu verfassen.
Es ist dir nicht erlaubt, auf Beiträge zu antworten.
Es ist dir nicht erlaubt, Anhänge hochzuladen.
Es ist dir nicht erlaubt, deine Beiträge zu bearbeiten.

BB-Code ist an.
Smileys sind an.
[IMG] Code ist an.
HTML-Code ist aus.
Trackbacks are an
Pingbacks are an
Refbacks are aus

Gehe zu:

Impressum · AGB · Datenschutz · Nach oben
Alle Zeitangaben in WEZ +1. Es ist jetzt 17:06 Uhr.
Powered by vBulletin® Copyright ©2000 - 2019, Jelsoft Enterprises Ltd.
LinkBacks Enabled by vBSEO © 2011, Crawlability, Inc.
Delphi-PRAXiS (c) 2002 - 2019 by Daniel R. Wolf