Имеется следующая задача по численным методам: найти собственное значение трехдиагональной матрицы А
методом парабол.
На пальцах алгоритм в общем то ясен, но с реализацией проблемы, ниже прикрепляю то что сделал сам, но в итоге
получается что после нескольких проходов на вход процедуры подаются три нуля, а этого быть не должно.
Просьба тем, кто хоть немного в этом разбирается, помочь найти и исправить ошибки.
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 невозможно, а значит невозможно найти и следующее приближение.