Решение СЛАУ методом LU-разложнеия Опции

[legend-muay]

Помогите найти ошибку в программе
При числах меньше 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.

Сообщение отредактировано: Lapp - 31.10.2009 9:23