Lineares Gleichungssystem lösen

[Kegasetu]

Ich habe mich mit den Operatoren von TMS etwas länger auseinandergesetzt, jedoch liegt der Entschluss, dass dies nicht unbedingt notwendig ist. In Excel ist es relativ einfach solche Gleichungssysteme zu lösen, dafür gibt es fast schon komplett vorgefertigte Funktionen.
Während eine Internetsuche bin ich bezüglich Gauß und Delphi dann auf folgenden Code gestoßen:

Delphi-Quellcode:

Beispiel 15.2 Der allgemeine Gauß'sche Algorithmus.
Diesmal verwenden wir ein dynamische Array. Das zweidimensionale Array wird deklariert durch den Typ:

type TarrayOfArrayOfExtended = array of array of extended;

Falls koeff[i,i] = 0 ist müssen noch die Spalten vertauscht werden. Das geschieht mit der Permutation p und Umkehrpermutation q.
Hier musst Du noch die Untit MathMohr miteinbinden. (oder die Entsprechenden Prozeduren durch eigene ersetzten.)

function IsInteger(const x: extended; eps: extended): boolean;
begin //eps globale Variable. Zum Beispiel eps = 1E-9
  result := frac(abs(x) + eps) < eps*2;
end;

FUNCTION ggtInt(a,b:longint):longint;
  begin if b=0 then result:=a
        else result:=ggtInt(b,a mod b);
  end;

FUNCTION ggtReal(a,b:Extended):Extended;
  begin if (a < maxlongint) and (b < maxlongint) then
           result:=ggTInt(round(a),round(b))
        else Begin
          if abs(b) < 0.5 then result:=a
          else result:=ggtReal(b,a-b*int(a/b));
        end;
end;

FUNCTION kgVReal(a, b: extended): extended;
begin result := a * b / ggTReal(a, b) end;

procedure LSG(n: integer; var aa: TarrayOfArrayOfExtended; var xx: array of extended);
             //Arrray aa[0..n,0..n+1]
var i0, j0,j0max: integer; //n Unbekannte bzw. Gleichungen
  p, q: array of integer; //permutation q=inv_p
  procedure invers; //erzeugt q = inverse Permutation von p
  var u, v: integer;
  begin
    for u := 0 to n - 1 do
      for v := 0 to n - 1 do
        if p[u] = v then q[v] := u;
  end;
  procedure tausche_sp(i, j: integer); // Spalten werden ausgetauscht
  var u, k: integer;
    x: extended;
  begin
    for u := 0 to n - 1 do Begin
      x := aa[u, i];
      aa[u, i] := aa[u, j];
      aa[u, j] := x; //=altes aa[u,i]
    End;
    k := p[i];
    p[i] := p[j];
    p[j] := k; //altes p[i]
    invers;
  end;

  procedure macheZeileGanzzahlig(zeile: integer); //bis auf rechte Seite  aa[z,n+1]
  var k : integer;
    d, zae, nen: extended; //wird das kgV des Nenners
  begin
    try
      d := 1;
      for k := 0 to n - 1 do Begin
        if not ErmittleBruch(abs(aa[zeile, k]), zae, nen,g_eps) then exit;
        d := kgVReal(nen, d);
      End;
      for k := 0 to n do aa[zeile, k] := d * aa[zeile, k];
      //Jetzt noch kürzen
      if not isInteger(aa[zeile, 1], g_eps) then exit;
      d := round(aa[zeile, 1]);
      for k := 0 to n - 1 do Begin
        if not isInteger(aa[zeile, k], g_eps) then exit;
        if d = 0 then d := round(aa[zeile, k]); //falls aa[zeile,1..]=0
        if abs(aa[zeile, k]) > 0 then
          d := ggTReal(round(aa[zeile, k]), d);
      End;
      if d <> 0 then for k := 0 to n do aa[zeile, k] := aa[zeile, k] / d;
    except {dann halt nicht} end;
    if aa[zeile, zeile] < 0 then for k := 0 to n do aa[zeile, k] := -aa[zeile, k];
  end;

  procedure VereinfacheRest(i: integer);
  var zeile, spalte: integer;
    d: extended;
  begin
    for zeile := 0 to n - 1 do if zeile <> i then Begin
        d := aa[zeile, i] / aa[i, i];
        if d <> 0 then Begin
          for spalte := 0 to n do if spalte <> i then
              aa[zeile, spalte] := aa[zeile, spalte] - d * aa[i, spalte] else
              aa[zeile, i] := 0; //=aa[zeile,i]-aa[zeile,i]/aa[i,i]*aa[i,i]
          macheZeileGanzzahlig(zeile);
        End;
      End;
  end;

begin //Hauptprogramm
  setlength(p,n+1);
  setlength(q,n);
  for j0 := 0 to n - 1 do Begin
    p[j0] := j0;
    q[j0] := j0;
  End;
  for i0 := 0 to n - 1 do Begin
    j0max := i0;
    for j0 := i0 + 1 to n - 1 do if abs(aa[i0, j0]) > abs(aa[i0, j0max]) then j0max := j0;
    if aa[i0,j0max] = 0 then Begin
      showmessage('Keine eindeutige Lösng!');
      exit;
    End;
    VereinfacheRest(i0);
  end;
  for i0 := 0 to n - 1 do xx[p[i0]] := aa[i0, n] ;
end;

procedure ArrayInMemo(n: integer; aa: TarrayOfArrayOfExtended; m:Tmemo);
  var i, j: integer;   //aa n zeilen und n + 1 Spalten
      s: string;
begin
  m.Lines.Clear;
  for i := 0 to n-1 do Begin
    s := '';
    for j := 0 to n do
      s := s + ' ' + floatToStr(aa[i,j]);
    m.Lines.Add(s);
  End;
end;

procedure TForm1.Button1Click(Sender: TObject);
  var n, i: integer;
      koeff: TarrayOfArrayOfExtended;
      a: array of extended;
begin
  n := 3; // von 0 bis 2
  setlength(a, n);
  setlength(koeff, length(a),length(a) + 1);
  koeff[0,0] := 5;  koeff[0,1] := -7; koeff[0,2] := 1; koeff[0,3] := 9;
  koeff[1,0] := -3; koeff[1,1] := 2;  koeff[1,2] := 3; koeff[1,3] := 4;
  koeff[2,0] := 2;  koeff[2,1] := 3;  koeff[2,2] := 4; koeff[2,3] := 0;
  ArrayInMemo(length(a),Koeff,memo1);
  LSG(length(a),koeff,a);
  memo1.lines.add('Lösung');
  for i := 0 to n - 1 do
    memo1.lines.add(ReellZuBruch(a[i]));
end;

(https://kilchb.de/lektionen1ff.php)

Ich bekomme das Kotzen so kompliziert kann das doch nicht sein?

Gegeben ist eine 10x10 Matrix, welche gelöst werden muss.