Inverse einer matrix

[swarley]

okay, ich war etwas übereilig. ihr könnt ja gar nicht wissen, was ich da noch alles fabriziert habe
Aufgabe ist es ein Unit namens "MyUnitMatrix" zu schreiben, in der mittels Unterprogrammen 2 Matrizen zu addieren, subtrahieren und multiplizieren sind. Zusätlich soll die Inverse ermittelt werden, und am besten das ganze auch nochmal überprüft werden (inverse mit ausgangsmatrix multiplizieren soll einheitsmatrix ergeben sozusagen). Soweit sogut.
ich hatte ein "urfassung" aus der entstammt auch der quelltext. Leider hat meiner Professorin das ganze optisch nicht gepasst (mal abgesehen davon, dass es nicht richtig funktioniert hat[was an dieser stelle heißt: programm wurde ausgeführt, inverse lieferte allerdings nicht das richtige ergebnis, die erste spalte war richtig, die restlichen enthielten werte, nicht die richtigen])
In dieser "urfassung" muss der Benutzer zu erst den rang seiner beiden Matrizen eingeben. anschließend mit dem button weiter erscheinen die Stringgrids der entsprechenden Größe, die angegeben wurde. Anschließend werden diese mit Zahlen gefüllt und je nach Button ("+", "-", "*") wurde überprüft ob die matrizen den bedingungen entsprechen, ansonsten müssen werte geändert werden. Unter einem der Stringgrids befand sich dann der button "inverse" dazu nun der quelltext aus der "hauptunit" sozusagen

Delphi-Quellcode:

procedure TForm2.ButtonInversClick(Sender: TObject);
var
  i, j, zeilen, spalten: Integer;
begin
  zeilen := StrToInt(Zeilenedit1.Text);
  spalten := StrToInt(Spaltenedit1.Text);
  if zeilen = spalten then
  begin
    for i := 0 to zeilen - 1 do
      for j := 0 to spalten - 1 do
      begin
        mat1[i, j] := StrToFloat(StringGridMatrix1.Cells[j, i]);
      end;
  end
  else
  begin
    ShowMessage('Zeilen-und Spaltenanzahl muss gleich groß sein!');
    exit;
  end;

  setlength(mathilf, zeilen);
  for i := 0 to zeilen - 1 do
  begin
    setlength(mathilf[i], 2 * spalten);
  end;

  invers(zeilen); // da spalten=zeilen reicht es aus, eins davon zu übergeben

  setlength(matinv, zeilen);
  for i := 0 to zeilen - 1 do
  begin
    setlength(matinv[i], zeilen);
  end;

  for i := 0 to zeilen - 1 do
  begin
    for j := 0 to zeilen - 1 do
    begin
      matinv[i, j] := mathilf[i, j + zeilen];
    end;
  end;

  StringGridMatrix3.Visible := true;
  labelergebnis.Visible := true;

  for i := 0 to zeilen - 1 do
  begin
    for j := 0 to zeilen - 1 do
    begin
      StringGridMatrix3.Cells[j, i] := FloatToStr(matinv[i, j]);
    end;
  end;
end;

end;

so, die matrizen, die eingegeben wurden also "mat1" und "mat2" habe ich global deklariert (was leider aufgabentechnisch auch falsch war). Ich hätte es lokal deklarieren und in das unterprogramm mit übergeben sollen. Das muss ich auch erneut überarbeiten.
naja und dann in meiner Unit "myunitmatrix" sozusagen, habe ich das unterprogramm deklariert, das is sozusagen der Quelltext oben. ich poste ihn trotzdem nochmal.

Delphi-Quellcode:

procedure invers(zeilen: Integer);
var
  i, j, k, pivotzeile: Integer;
  mateinheit: array of array of double;
  pivotelement, faktor, hilf: double;
begin
  setlength(mateinheit, zeilen);
  for i := 0 to zeilen - 1 do
  begin
    setlength(mateinheit[i], zeilen);
    mateinheit[i, i] := 1;
  end;

  for i := 0 to zeilen - 1 do
  begin
    for j := 0 to zeilen - 1 do
    begin
      mathilf[i, j] := mat1[i, j];
    end;
  end;

  for i := 0 to zeilen - 1 do
  begin
    for j := 0 to zeilen - 1 do
    begin
      mathilf[i, j + zeilen] := mateinheit[i, j];
    end;
  end;

  for k := 0 to zeilen - 2 do
  begin
    pivotelement := 0;
    pivotzeile := k;
    for i := k to zeilen - 1 do
    begin
      if abs(mathilf[i, k]) > pivotelement then
      begin
        pivotelement := abs(mathilf[i, k]);
        pivotzeile := i;
      end;
      if k <> pivotzeile then
      begin
        for j := k to zeilen do
        begin
          hilf := mathilf[k, j];
          mathilf[k, j] := mathilf[i, j];
          mathilf[i, j] := hilf;
        end;
      end;
    end;
    if mathilf[k, k] = 0 then
      exit;
    for i := k + 1 to zeilen - 1 do
    begin
      faktor := mathilf[i, k] / mathilf[k, k];
      for j := k to zeilen do
      begin
        mathilf[i, j] := mathilf[i, j] - faktor * mathilf[k, j];
      end;
    end;
  end;

  for k := zeilen - 1 downto 0 do
  begin
    if mathilf[k, k] = 0 then
      exit;
    for i := k - 1 downto 0 do
    begin
      faktor := mathilf[i, k] / mathilf[k, k];
      for j := zeilen downto i do
      begin
        mathilf[i, j] := mathilf[i, j] - faktor * mathilf[k, j];
      end;
    end;
    for j := zeilen to 2 * zeilen - 1 do
    begin
      mathilf[k, j] := mathilf[k, j] / mathilf[k, k];
    end;
  end;
end;

ich weiß dieser beitrag ist mittlerweile schon überlang, aber ich komm echt nicht mehr klar. ich bräuchte irgendwie einen funktionierenden algorithmus für diese gauß jordan geschichte. Ich habe dummerweiße (glaub ich) auch noch nicht den fall abgegriffen, was geschieht wenn das erste element der ersten zeile eine 0 ist.
ich wäre auch über neue denkanstöße oder hilfelinks sehr dankbar, denn da ich selber nicht mehr checke was ich da hin implementiert habe, besstehe ich nicht auf das gedöns ^^
vielen dank an alle die sich die zeit nehmen das anzusehen

PS: ich habe die beschriebenen dateien mal als zip angehängt falls ich in meiner ausführung zu ungenau oder schlusig war.

[BUG]

Zitat von swarley:

ich wäre auch über neue denkanstöße oder hilfelinks sehr dankbar, denn da ich selber nicht mehr checke was ich da hin implementiert habe, besstehe ich nicht auf das gedöns ^^

Wie wäre es, wenn du die Zeilenoperationen in eigene Funktionen schiebst?
Das würde deinen Code viel lesbarer machen.


Du kannst Gauß "rekursiv" anwenden, wobei du immer eine Spalte und eine Zeile weniger betrachtest:

Angenommen du hast eine n*n-Matrix M = (a_ij) mit i = Zeilennummer, j = Spaltennummer.

Code:

für j := 1 bis n tue
  wenn a_jj = 0 dann
    suche i > j mit a_ij <> 0; // existiert, da Matrix vollen Rang hat
    tauche a_1j und a_ij;
  // jetzt ist a_jj <> 1
  Multipliziere Zeile j mit dem multiplikativen Inversen von a_jj;
  // jetzt ist a_jj = 1
  Addiere Zeile j so auf die restlichen Zeilen, dass dort Nullen in Spalte j entstehen (passenden Faktor wählen);

XXXX 1. Schleifendurchlauf
XXXX 2. Schleifendurchlauf
XXXX 3. Schleifendurchlauf
XXXX 4. Schleifendurchlauf

Apropos, in welchem Körper rechnest du eigentlich?

[swarley]

stimmt, an rekursion habe ich noch gar nicht gedacht... werde mich daran mal versuchen.
also in der aufgabenstellung ist nichts gegeben, demnach nehm ich ma ganz stark an, dass wir uns in Q befinden.
es soll grob gesehn ein Programm für Schüler entstehen, demnach nehme ich an, dass eben einfach positive und negative Dezimalzahlen mit und ohne Nachkommastellen eingegeben werden sollen
wie schaffe ich es denn, dass ich die eingelesen matrizen dem unterprogramm übergeben kann. muss ich dazu eine typdeklaration machen? also ich meine:
type Matrizen = TMatrix usw? wenn ich die arrays als lokale variablen einfach übergeben, wird mir das rot unterwellt.

[BUG]

Zitat von swarley:

also in der aufgabenstellung ist nichts gegeben, demnach nehm ich ma ganz stark an, dass wir uns in Q befinden.

Das solltest du unbedingt klären oder den Algorithmus unabhängig vom Körper implementieren (aber deinem bisher gezeigtem Quellcode nach solltest du das besser nicht versuchen).

Die reellen Zahlen solltest du vermeiden, da man die im Computer nicht genau darstellen kann und deshalb nur Ärger hat.
Wenn du es dir aussuchen darfst: mit einer Restklasse modulo einer Primzahl gibt es imho am wenigsten Probleme.

[Bjoerk]

Gegen eine Typendeklaration spricht in Delphi eigentlich nie etwas.

Delphi-Quellcode:

type
  TMat = array of array of double;
  TVek = array of double;

Vielleicht noch ein kleiner Hinweis. Die inverse Matrix zu berechnen ist sehr anspruchsvoll. Deshalb geht man in der Regel anders vor: N mal Gauss mit Einheitsvektor als Rechter Seite durchlaufen ergibt (Spaltenweise) die Reziproke Matrix.

Beispiel:

Matrix:
1 2 0
2 3 0
3 4 1

Rechte Seite 1. Durchlauf:
1
0
0
Lösungsvektor 1. Durchlauf:
-3
2
1

Rechte Seite 2. Durchlauf:
0
1
0
Lösungsvektor 2. Durchlauf:
2
-1
-2

Rechte Seite 3. Durchlauf:
0
0
1
Lösungsvektor 3. Durchlauf:
0
0
1

Reziproke Matrix:
-3 2 0
2 -1 0
1 -2 1

[BUG]

Zitat von Bjoerk:

Vielleicht noch ein kleiner Hinweis. Die inverse Matrix zu berechnen ist sehr anspruchsvoll. Deshalb geht man in der Regel anders vor: N mal Gauss mit Einheitsvektor als Rechter Seite durchlaufen ergibt (Spaltenweise) die Reziproke Matrix.

Lustigerweise entspricht das genau dem Algorithmus, über den wir die ganze Zeit reden (nur das man da alle Durchläufe gleichzeitig macht)

[Bjoerk]

Hi Bug,

stimmt, auf die Idee bin ich noch gar nicht gekommen. Man kann den Gauss ja logischerweise für mehrere Reche Seiten gleichzeitig laufen lassen. Danke für den Hinweis.

Gruß
Thomas