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