При числах меньше 10 работает нормально, если есть числа больше 10 начинает работать неправильно
matr=array[1..100,1..100] of extended;
TResult=array[1..100]of extended;
const
eps=1e-12;
{
Раскаладывает матрицу А на произведение LU
А= L11 U12 U13 {n=3}
L21 L22 U23
L31 L32 L33
}
function LUDecomposition(var a:matr;n:word):boolean;
var i,j,k:word;
s:extended;
u,l:matr;
begin
if abs(a[1,1])<eps then
begin
LUDecomposition:=false;
exit;
end;
for i:=2 to n do
for k:=2 to n do
begin
if abs(a[i,i])<eps then
begin
LUDecomposition:=false;
exit;
end;
a[1,k]:=a[1,k]/a[1,1];
if k<=i then
begin
s:=0;
for j:=1 to k-1 do
s:=s+a[i,j]*a[j,k];
a[i,k]:=a[i,k]-s;
end
else
begin
s:=0;
for j:=1 to i-1 do
s:=s+a[i,j]*a[j,k];
a[i,k]:=(a[i,k]-s)/a[i,i];
end;
end;
if abs(a[n,n])<eps then
begin
LUDecomposition:=false;
exit;
end;
LUDecomposition:=true;
end;
{
раскладываем A на L и U
L= L11 0 0
L21 L22 0
L31 L32 L33
U= 1 U12 U13
0 1 U23
0 0 1
}
procedure Decomposition(a:matr;n:byte;var NewL,NewU:matr);
var i,j:word;
begin
for i:=1 to n do
for j:=1 to n do
begin
if i>=j
then
begin
NewL[i,j]:=a[i,j];
if i=j then NewU[i,j]:=1
else NewU[i,j]:=0;
end
else
begin
NewU[i,j]:=a[i,j];
NewL[i,j]:=0;
end;
end;
end;
{
Обратный ход для
Ly=b
Ux=y
}
procedure BackStep(a:matr;b:TResult;n:byte;flag:char;var res:TResult);
var i,j,k:word;
s:extended;
begin
if flag in['U','u'] then
begin
res[n]:=b[n]/a[n,n];
for i:=n-1 downto 1 do
begin
s:=0;
for j:=i+1 to n do
s:=s+(a[i,j]*res[j]);
res[i]:=(b[i]-s)/a[i,i];
end;
end;
if flag in['L','l'] then
begin
res[1]:=b[1]/a[1,1];
for i:=2 to n do
begin
s:=0;
for j:=1 to i-1 do
s:=s+(a[i,j]*res[j]);
res[i]:=(b[i]-s)/a[i,i];
end;
end;
end;
{
Основная ПП
}
begin
read_matrix(n,a,b);{считываем матрицу и столбец свободных членов}
if LUDecomposition(a,n)=false then exit;
Decomposition(a,n,l,u);
BackStep(l,b,n,'l',y);
BackStep(u,y,n,'u',x);
end.