Метод парабол поиска собственных значений матрицы Опции

[redeezko]

Всем доброго времени суток!
Имеется следующая задача по численным методам: найти собственное значение трехдиагональной матрицы А
методом парабол.

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

program Project2;
{$APPTYPE CONSOLE}
uses
  SysUtils;
const n=3; e=0.01;
type
 tmatrix=array[1..n,1..n] of real;
 tvector=array[0..n] of real;
var l,i,j:integer;
l_0,l_1,l_2,det1,det2,det3,xn1,xn:real;
a:tmatrix; d:tvector;

procedure matr(var a:tmatrix);
var i:integer;
begin
  for i:=1 to n do
  a[i,i]:=random(10-1)+1;
 for i:=1 to n-1 do
  a[i,i+1]:=random(10-1)+1;
 for i:=2 to n do
  a[i,i-1]:=random(10-1)+1;
end;

function det(l:real; a:tmatrix):real;
var i,j:integer;

begin
 for i:=1 to n do
  a[i,i]:=a[i,i]-l;
 d[0]:=1;
 d[1]:=a[1,1];
 for i:=2 to n do
  d[i]:=(a[i,i]*d[i-1])-(a[i,i-1]*a[i-1,i]*d[i-2]);
 det:=d[n];
end;

procedure parabol(var x1,x2,x3,f1,f2,f3:real);
var a1,b1,c1,x11,x33,xp1,xp2,z:real;
begin
x11:=x1-x2;
x33:=x3-x2;
if (x11=0) or (x11=x33) then
begin
 a1:=0;
 b1:=(f3-f2)/x33;
 c1:=f2;
 xp1:=-c1/x33;
 xp2:=xp1;
end;
if (x33=0) or (x11=x33) then
begin
 a1:=0;
 b1:=(f1-f2)/x11;
 c1:=f2;
 xp1:=-c1/x11;
 xp2:=xp1;
end;
if (x11<>0) and (x33<>0) and (x11<>x33) then
begin
z:=x11*x11*x33-x33*x33*x11;
a1:=((f1-f2)*x33-(f3-f2)*x11)/z;
b1:=(f1-f2-a1*x11*x11)/x11;
c1:=f2;
end;
if a1<>0 then begin
xp1:=((-1)*b1-sqrt(b1*b1-4*a1*c1))/(2*a1);
xp2:=((-1)*b1+sqrt(b1*b1-4*a1*c1))/(2*a1)
end;

if abs(xp1-x3)<abs(xp2-x3) then xn:=xp1
 else
  xn:=xp2;
if abs(xn-xn1)>e then
begin
xn1:=xn;
if xn=xp1 then
 begin
  f3:=det(xn,A);
  parabol(xn,x1,x2,f3,f1,f2);
  end
 else
   begin
    f1:=det(xn,A);
    parabol(x2,x3,xn,f2,f3,f1);
   end;
end;
end;

begin
randomize;
matr(a);
 for i:=1 to n do
  begin
   writeln;
   for j:=1 to n do
    write(a[i,j]2,' ');
  end;
l_0:=-1;
l_1:=0;
l_2:=1;
det1:=det(l_0,a);
det2:=det(l_1,a);
det3:=det(l_2,a);
xn1:=0;
writeln;
parabol(l_0,l_1,l_2,det1,det2,det3);
writeln(xn1:4:2);
readln;

end.

как видно из кода, основные входные значения - матрица, задается рандомно (ибо метод должен работать для любой матрицы)
l0,l1,l2 - исходные 3 точки: -1, 0 ,1 - стандартно для метода парабол.
Проблема - после нахождения энного (каждый раз разного - зависит от матрицы) приближения случается такая ситуация что х11,х33 становятся одновременно равными нулю, а соответственно решить систему относительно a1,b1 невозможно, а значит невозможно найти и следующее приближение.



Сообщение отредактировано: redeezko - 21.12.2011 20:01