Lineares Gleichungssystem lösen?

[paresy]

ich habe mir den algo aus der code library angeschaut und musste leider feststellen, dass er irgendwie bei mir nicht funktioniert :(

beispiel:

Delphi-Quellcode:

program gauss;

{$APPTYPE CONSOLE}

uses
  SysUtils, Math;

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;
    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, 4, 5);
  A[0][0] := 0;  A[0][1] := 0;  A[0][2] := 0;   A[0][3] := 1;  A[0][4] := 0;
  A[1][0] := 1;  A[1][1] := 1;  A[1][2] := 1;   A[1][3] := 1;  A[1][4] := 1;
  A[2][0] := 8;  A[2][1] := 4;  A[2][2] := 2;   A[2][3] := 1;  A[2][4] := 4;
  A[3][0] := 64; A[3][1] := 16; A[3][2] := 4;   A[3][3] := 1;  A[3][4] := 8;

  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, 4, 5);
  WriteLn;

  for i := 0 to High(Res) do
    WriteLn('x(', i + 1, ') = ', FloatToStr(Res[i]));

  ReadLn;
end.

die lösung die ich auf dem papier errechnet habe ist: -1/3*X^3+2*X^2-2/3*X

ich glaube, dass es am gauss algoritmus liegt, weil es in mehreren anderen programmen die den algo nutzen nicht geht.

hat einer vielleicht einen alternativen algo?

grüße, paresy

[Jelly]

Zitat:

Zitat von paresy

die lösung die ich auf dem papier errechnet habe ist: -1/3*X^3+2*X^2-2/3*X

1) Wie sieht dein Ausgangsproblem aus, dein lineares Gleichungssystem
2) Wie soll ich deine genannte Lösung interpretieren... Da seh ich lediglich einen math. Ausdruck
3) Was liefert dir denn der Algo hier aus dem Forum
4) Bist du sicher, daß du ein lineares Gleichungssystem hast... Bei dir tauchen Potenzen auf.

Ohne die Informationen kann hier keiner was zu deinem Problem sagen.

[Dust Signs]

Man sieht auf den ersten Blick, dass dein Gleichungssystem nicht linear ist. Gauss bzw. Gauss-Jordan können aber nur mit linearen Gleichungssystemen rechnen. Außerdem: darf ich fragen, warum deine Lösung wie eine kubische Gleichung aussieht?

btw: hier mein Gauss-Jordan, is etwas effektiver als Gauss:

//EDIT: musst den Code natürlich noch anpassen; is ein älterer Code, an dem ich nicht mehr arbeite

Delphi-Quellcode:

program gauss_jordan;

{$APPTYPE CONSOLE}

uses
  SysUtils;

var
  i, j, k: Integer; //für Schleifen
  pivot, z: Integer; //für Pivotelement und gemerktes Element; zum Multiplizieren

var
  matrix: Array [1..3] of Array [1..3] of Integer;
  loesungsvektor: Array [1..3] of Integer;
  ergebnis: Array [1..3] of Integer;

procedure ZwischenergebnisAnzeigen;
var
  l, m: Integer;
begin
  WriteLn;
  for l := 1 to 3 do begin
    WriteLn;
    for m := 1 to 3 do begin
      Write(matrix[l, m]: 12);
      end;
    Write('  |', loesungsvektor[l]: 12);
    end;
  WriteLn;
  WriteLn('z: ', z, ' pivot: ', pivot, ', j: ', j, ', k: ', k);
end;

begin

  WriteLn('Gleichungslösen nach Gauss-Jordan');
  WriteLn('=================================');
  WriteLn;

  WriteLn('Bitte geben Sie die Komponenten der Koeffizientenmatrix ein (Zeile, Spalte):');
  WriteLn;

  for i := 1 to 3 do begin
    for j := 1 to 3 do begin
      Write('Element ', i, ',', j, ': ');
      ReadLn(matrix[i, j]);
      end;
    end;

  WriteLn;
  WriteLn('Bitte geben Sie die Komponenten des Lösungsvektors ein:');
  WriteLn;

  for i := 1 to 3 do begin
    Write('Komponente ', i, ': ');
    ReadLn(loesungsvektor[i]);
    end;

  //Gauss-Jordan

  for i := 1 to 3 do begin
    pivot := matrix[i, i]; //Pivotelement
    for j := 1 to 3 do begin //Zeilen
      if i <> j then begin //Pivotzeile unverändert lassen
        z := matrix[j, i]; //aktuelles Element merken
        for k := 1 to 3 do begin //Spalten
          ZwischenErgebnisAnzeigen;
          ReadLn;
          matrix[j, k] := pivot * matrix[j, k] - z * matrix[i, k]; //Elemente berechnen
          end;
        loesungsvektor[j] := pivot * loesungsvektor[j] - z * loesungsvektor[i]; //Komponenten des Lösungsvektors analog berechnen
        end;
      end;
    end;

    ZwischenergebnisAnzeigen;

  //Gauss-Jordan Ende

  for i := 1 to 3 do
    ergebnis[i] := loesungsvektor[i] div matrix[i][i];

  WriteLn;
  WriteLn;
  WriteLn('Lösung der Gleichung:');
  WriteLn;

  for i := 1 to 3 do
    WriteLn('Unbekannte ', i, ': ', ergebnis[i]);

  WriteLn;
  WriteLn('-----------------------------------------------------------------------------');
  WriteLn('Gleichungslösen nach Gauss-Jordan - (c) by Dust Signs Andreas Unterweger 2004');
  ReadLn;

end.

Dust Signs

//EDIT2+3: immer dieser Tippfehler ^^

[axelf98]

Ich hab mich auch schon mal mit der Materie beschäftigt...

Gleichungssysteme lösen

[paresy]

Zitat:

Zitat von axelf98

Ich hab mich auch schon mal mit der Materie beschäftigt...

Gleichungssysteme lösen

ja genau.

ich wollte halt aus 4 punkten eine gleichung ermitteln.

und dann, wie man es in mathe gelernt hat einsetzen und dann hat man 4 gleichungen.

diese dachte, ich man könnte die einfach so einsetzten und fertig wär die sache :)

quasi:

P(0/0)
P(1/1)
P(2/4)
P(4/8)

dann:

Code:

0 = 0  0  0  1
1 = 1  1  1  1
4 = 2³ 2² 2  1
8 = 4³ 4² 4  1

das in die matrix stecken und mit den algo ausrechnen lassen, sodass rauskommt

Code:

a = 1/3
b = 2
c = -2/3
d = 0

welches eingesetzt die oben genannte funktion ergäbe.

grüße, paresy

[axelf98]

Zitat:

Zitat von paresy

man könnte die einfach so einsetzten und fertig wär die sache :)

Theoretisch geht das auch, aber der Gaußalgorithmus mag ohne Modifikation nicht jedes Gleichungssystem.. Das war auch der Grund, warum ich mir meinen eigenen Algorithmus geschrieben habe :)

[paresy]

Zitat:

Zitat von axelf98

Zitat:

Zitat von paresy

man könnte die einfach so einsetzten und fertig wär die sache :)

Theoretisch geht das auch, aber der Gaußalgorithmus mag ohne Modifikation nicht jedes Gleichungssystem.. Das war auch der Grund, warum ich mir meinen eigenen Algorithmus geschrieben habe :)

so hab mir mal deine klasse angeguckt, jedoch weiß ich nicht was ich als l/r parameter übergeben soll. ( funktion: Systemloesen )

könntest du mir vielleicht nen tipp geben ;) ? bzw vllt eine kleine demo?

grüße, paresy

[axelf98]

Zitat:

Zitat von paresy

könntest du mir vielleicht nen tipp geben ;) ? bzw vllt eine kleine demo?

Kein Problem.. Ich hab aber gesehen, dass diese Parameter gar nicht gebraucht werden.. Wohl aber die Backup-Funktion, die man vorher füllen muss, damit am Ende das System noch mal geprüft werden kann... Unten mal ein Anwendungsbeispiel mit Editfeldern...

Die Unit ist nebenbei gesagt noch nicht perfekt, aber löst so gut wie alles, was lösbar ist und hat keine Probleme wie Gauß!

Delphi-Quellcode:

var LoesungSys : Loesungssystem;
    DieMatrix : TMatrix;
    Loesungen : TLoesung;
....

  setlength(DieMatrix,0,0);
  setlength(DieMatrix,Gleichungen,Unbekannte+2);  // vorbereiten

 // Daten aus Edits holen...
 for i := 0 to Gleichungen-1 do
  for j := 0 to Unbekannte do
    begin
    try
     if Felder[i,j].Text <> '' then
     DieMatrix[i,j] := StrtoFloat(Felder[i,j].Text) ELSE
     DieMatrix[i,j] := 0;
     if j = unbekannte then
     begin
      DieMatrix[i,j] :=  - DieMatrix[i,j];
      if DieMatrix[i,j] = 0 then
      begin
       DieMatrix[i,j] := 1E-1000;  // Rundungsfehler wegen der Null
       genau := false;
      end;
     end;
    except
     MessageDLG('Fehler bei der Eingabe in Zeile '+inttostr(i+1)+', Spalte '+ inttostr(j+1),mtinformation,[mbOK],0);
      Felder[i,j].SetFocus;
     Felder[i,j].SelectAll;
     exit;
    end;
   end;
   self.LoesungSys.backup(DieMatrix,Gleichungen,unbekannte+1);  // WICHTIG ZUM ÜBERPRÜFEN SPÄTER
   setlength(loesungen,unbekannte);
                // Lösen...
   try
  Loesungen := LoesungSys.Systemloesen(DieMatrix,Gleichungen,unbekannte+1,unbekannte,1,true);       // Lösen...
   except
   MessageDLG('Das Gleichungssystem ist nicht lösbar!',mtinformation,[mbOK],0);
   self.StatusBar.Panels[1].Text := 'System nicht lösbar.';
   abort;
   end;

   self.meldungen := self.LoesungSys.Meldungen; // Fehler und Meldungen

[paresy]

Zitat:

DieMatrix[i,j] := - DieMatrix[i,j];
if DieMatrix[i,j] = 0 then
begin
DieMatrix[i,j] := 1E-1000; // Rundungsfehler wegen der Null
genau := false;
end;

das war ne spannende konstruktion :)

hat gut funktioniert, aber irgendwie bekam ich immer ne access violation wenn ich das programm geschlossen habe, welche ich einfach nicht ausfindig machen konnte.

hab aber noch etwas anderes gefunden, was auch recht gut zu sein scheint:

http://www.gkinf.de/sdhp/linglsys.html

grüße, paresy

[axelf98]

Zitat:

Zitat von paresy

hat gut funktioniert, aber irgendwie bekam ich immer ne access violation wenn ich das programm geschlossen habe, welche ich einfach nicht ausfindig machen konnte.

Ja, etwas abendteuerlich ist es schon, aber es läuft... Das mit der Violation liegt wahrscheinlich daran, dass du irgendeinen Speicherbereich zu klein initialisiert hast. Das Problem hatte ich auch gehabt. Also mal alle Setlengths untersuchen...

Gruß