LU-разложение (Delphi)

[marwell.]

доброго времени суток
цель: реализовать метод LU-разложения для решения СЛАУ размерности 6х6
метод несложный, вроде все верно делаю, а корни неправильные получаются
сама матрица передается в ish_matr, свободные члены в svob_koeffs

  type
    Tmass=array[1..100] of double;
    Tmatrix=array[1..100,1..100] of double;
var
  Form1: TForm1;
  korni,z:Tmass;
  ish_matr:TMatrix; svob_koeffs:TMass;
      L,U: TMatrix;

const
krat=6;

implementation

{$R *.dfm}

procedure TForm1.Button1Click(Sender: TObject);

procedure LU(a:TMatrix; free:TMass);

function sum_korni(i,f:integer):double;
var j:integer;
    s:double;
begin
s:=0;
case f of
0: for j:=1 to i-1 do s:=s+z[j]*L[i,j];
1: for j:=i+1 to krat do s:=s+korni[j]*U[i,j];
end;
Result:=s;
end;

var
    i,j,k:integer;
begin
for i:=1 to krat do
for j:=1 to krat do U[i,j]:=a[i,j];
for i:=1 to krat do begin
for j:=i to krat do begin
L[j,i]:=U[j,i]/U[i,i];
end;
end;
for k:=2 to krat do begin
for i:=k-1 to krat do
for j:=i to krat do
L[j,i]:=U[j,i]/U[i,i];
for i:=k to krat do
for j:=k-1 to krat do
U[i,j]:=U[i,j]-L[i,k-1]*U[k-1,j];
end;
z[1]:=free[1]/L[1,1];
for i:=2 to krat do begin
z[i]:=(free[i]-sum_korni(i,0))/L[i,i];
end;
korni[krat]:=z[krat];
for i:=krat-1 downto 1 do begin
korni[i]:=z[i]-sum_korni(i,1);
end;
end;


var i, j: integer;
begin
...
  LU(ish_matr,svob_koeffs);
  Memo1.Clear;
    for i:=1 to 6 do Memo1.Lines.Add('p['+IntToStr(i)+']='+FloatToStrf(korni[i], ffGeneral, 8, 3));
end;

end.

P.S сами матрицы L и U проверял, верно составляются

[Sciv]

for i:=1 to krat do
for j:=1 to krat do U[i,j]:=a[i,j];

Я думаю, вполне прокатит u:=a;

L[j,i]:=U[j,i]/U[i,i];

Странно, что правильно формирует L и U. U[j,i] у Вас нигде до этой операции не формируется, а по идее должно быть равно:

for k:=1 to i-1 do sum:=sum+l[j,k]*u[k,i];
u[j,i]:=a[j,i]-sum;

P.S сами матрицы L и U проверял, верно составляются

И при перемножении дают исходную?
Прикрепите проект к посту, так проще будет разбираться, код у Вас трудновато читается.

[marwell.]

Странно, что правильно формирует L и U. U[j,i] у Вас нигде до этой операции не формируется, а по идее должно быть равно:

for i:=1 to krat do
for j:=1 to krat do U[i,j]:=a[i,j];

тут в массив u передается исходный массив a

архив прикрепил

[marwell.]

Странно, вроде все правильно делаю

[Sciv]

Ну вот что у меня из Вашей программы вышло (для проверки разложения я там много чего добавил).

Вывод: при первом нажатии на кнопку корни считаются правильно. Но вот дальше - это непонятная проблема, которую я списал на точность числа с плавающей точкой (они там получаются при разложении) - потому что ничего более адекватного и логичного не придумал. Судите сами: исходные данные - одни и те же, операции над ними тоже, а результат разный. Причем( если верить трассировке) до процедуры LU все нормально, а вот после - расхождения...

[marwell.]

спасибо за ответ! а вывод действительно странный