IPB
ЛогинПароль:

> Прочтите прежде чем задавать вопрос!

1. Заголовок темы должен быть информативным. В противном случае тема удаляется ...
2. Все тексты программ должны помещаться в теги [code=pas] ... [/code], либо быть опубликованы на нашем PasteBin в режиме вечного хранения.
3. Прежде чем задавать вопрос, см. "FAQ", если там не нашли ответа, воспользуйтесь ПОИСКОМ, возможно такую задачу уже решали!
4. Не предлагайте свои решения на других языках, кроме Паскаля (исключение - только с согласия модератора).
5. НЕ используйте форум для личного общения, все что не относится к обсуждению темы - на PM!
6. Одна тема - один вопрос (задача)
7. Проверяйте программы перед тем, как разместить их на форуме!!!
8. Спрашивайте и отвечайте четко и по существу!!!

 
 Ответить  Открыть новую тему 
> Метод парабол поиска собственных значений матрицы
сообщение
Сообщение #1


Новичок
*

Группа: Пользователи
Сообщений: 22
Пол: Мужской

Репутация: -  0  +


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

На пальцах алгоритм в общем то ясен, но с реализацией проблемы, ниже прикрепляю то что сделал сам, но в итоге
получается что после нескольких проходов на вход процедуры подаются три нуля, а этого быть не должно.
Просьба тем, кто хоть немного в этом разбирается, помочь найти и исправить ошибки.
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]4.gif2,' ');
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 -
 Оффлайн  Профиль  PM 
 К началу страницы 
+ Ответить 
сообщение
Сообщение #2


Гуру
*****

Группа: Пользователи
Сообщений: 1 013
Пол: Мужской
Ада: Разработчик
Embarcadero Delphi: Сторонник
Free Pascal: Разработчик

Репутация: -  627  +


Просьба к тем, кто выкладывает код - либо выкладывать его полностью, либо (хотя бы) говорить, с какими значениями входных переменных (ВСЕХ, ибо есть еще и глобальные, насколько я вижу) он тестировался, когда было замечено вышеописанное поведение.
 Оффлайн  Профиль  PM 
 К началу страницы 
+ Ответить 
сообщение
Сообщение #3


Новичок
*

Группа: Пользователи
Сообщений: 22
Пол: Мужской

Репутация: -  0  +


отредактировал. что еще непонятно могу пояснить
 Оффлайн  Профиль  PM 
 К началу страницы 
+ Ответить 
сообщение
Сообщение #4


Гуру
*****

Группа: Пользователи
Сообщений: 1 013
Пол: Мужской
Ада: Разработчик
Embarcadero Delphi: Сторонник
Free Pascal: Разработчик

Репутация: -  627  +


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


Эскизы прикрепленных изображений
Прикрепленное изображение
 Оффлайн  Профиль  PM 
 К началу страницы 
+ Ответить 
сообщение
Сообщение #5


Новичок
*

Группа: Пользователи
Сообщений: 22
Пол: Мужской

Репутация: -  0  +


Цитата
Не получается этого. Получается, что в процедуре происходит деление на 0, поэтому она завершается аварийно:

Да, немного неверно выразился. Подаются 3 одинаковых числа, и в итоге x11 и х33 становятся равными 0, а на них делить надо..

Сообщение отредактировано: redeezko -
 Оффлайн  Профиль  PM 
 К началу страницы 
+ Ответить 
сообщение
Сообщение #6


Новичок
*

Группа: Пользователи
Сообщений: 22
Пол: Мужской

Репутация: -  0  +


Нашел решение данной задачи на С++, но в нем практически не разбираюсь.
Буду очень признателен, если кто нибудь поможет перевести данный код на Паскаль.

//---------------------------------------------------------------------------

#pragma hdrstop

//---------------------------------------------------------------------------
#include<iostream.h>
#include <iomanip>
#include<math.h>
#include<windows.h>
#include <stdlib.h>
#pragma argsused
const double eps=1;
const int n=6;
double det(double *a, double *b, double *c, double x)
{
double b1[n];
for(int i=0;i<n;i++)
b1[i]=b[i];
b1[0]-=x;
double rez=b1[0];
for(int i=1;i<n;i++)
{
b1[i]-=x;
if(b1[i-1]==0){cout<<"ошибка: деление на 0 при вычислении определителя!!! строка "<<i<<" при x="<<x<<endl; system("pause"); exit(0);}
b1[i]+=-c[i-1]*a[i-1]/b1[i-1];
rez*=b1[i];
}
return rez;
}
int main(int argc, char* argv[])
{
SetConsoleCP(1251);
SetConsoleOutputCP(1251);
randomize();
double a[n-1], b[n], c[n-1];
double x1=-1, x2=0, x3=1;

//cout<<"введите первую диагональ: "<<endl;
for(int i=0;i<n-1;i++)
//cin>>a[i];
a[i]=random(20)-10;
//cout<<"введите вторую(главную) диагональ: "<<endl;
for(int i=0;i<n;i++)
//cin>>b[i];
b[i]=random(20)-10;
//cout<<"введите третью диагональ: "<<endl;
for(int i=0;i<n-1;i++)
//cin>>c[i];
c[i]=random(20)-10;

cout<<endl<<setprecision(4)<<b[0]<<" "<<c[0];
for(int i=0; i<n-2; i++) cout<<" 0";
cout<<endl;
for(int i=1;i<n-1;i++)
{
for(int j=1;j<i;j++)cout<<"0 ";
cout<<a[i-1]<<" "<<b[i]<<" "<<c[i]<<" ";
for(int j=0;j<n-i-2;j++)cout<<"0 ";
cout<<endl;
}
for(int j=1;j<n-1;j++)cout<<"0 ";
cout<<a[n-2]<<" "<<b[n-1]<<" ";
cout<<endl;
//-------------------- первое соственное значение --------
double p, q, f1, f2, f3, xRez1, xRez2, rez;

do
{
f1=det(a,b,c,x1);
f2=det(a,b,c,x2);
f3=det(a,b,c,x3);
if((x1-x2)*(x3-x2)*(x1-x3)==0)
{cout<<"ошибка: деление на 0!!!"; system("pause"); return 0;}
p=((x3-x2)*(f1-f2)-(x1-x2)*(f3-f2))/((x1-x2)*(x3-x2)*(x1-x3));
if((x1-x2)==0)
{cout<<"ошибка: деление на 0!!!"; system("pause"); return 0;}
q=(f1-f2-p*(x1-x2)*(x1-x2))/(x1-x2);
if(q*q-4*p*f2<0)
{
cout<<"метод сечений"<<endl;
while(abs(f3)>eps)
{
if(abs(f1-f3)<eps){cout<<"комплексные значения"<<endl;return 0;}
x3=-f1*(x1-x3)/(f1-f3)+x1;
f3=det(a,b,c,x3);
}
cout<<"первое собственное значение: "<<x3<<endl;
cout<<det(a,b,c,x3)<<endl;
rez=x3;
goto G;
}
{
xRez1=(-q+sqrt(q*q-4*p*f2))/(2*p);
xRez2=(-q-sqrt(q*q-4*p*f2))/(2*p);

if(abs(xRez2-x2)<abs(xRez1-x2)) xRez1=xRez2;
if(xRez1<x2)x1=xRez1;
else x3=xRez1;
}
}
while(abs(det(a,b,c,xRez1))>eps);
cout<<endl<<setprecision(36)<<"первое собственное значение: "<<xRez1<<endl;
cout<<det(a,b,c,xRez1)<<endl;
rez=xRez1;
//----------------- вычисляем второе собст значение
G:
x1=-1; x2=0; x3=1;
do
{
if(rez==x1||rez==x2||rez==x3)
{cout<<"деление на 0!!!"; system("pause");return 0;}
f1=(det(a,b,c,x1))/(x1-rez);
f2=det(a,b,c,x2)/(x2-rez);
f3=det(a,b,c,x3)/(x3-rez);
if((x1-x2)*(x3-x2)*(x1-x3)==0)
{cout<<"ошибка: деление на 0!!!"; system("pause"); return 0;}
p=((x3-x2)*(f1-f2)-(x1-x2)*(f3-f2))/((x1-x2)*(x3-x2)*(x1-x3));
if((x1-x2)==0)
{cout<<"ошибка: деление на 0!!!"; system("pause"); return 0;}
q=(f1-f2-p*(x1-x2)*(x1-x2))/(x1-x2);
if(q*q-4*p*f2<0)
{
cout<<"метод сечений"<<endl;
while(abs(f2)>eps)
{
if(abs(f1-f3)<eps){cout<<"комплексные значения"<<endl;return 0;}
x2=-f1*(x1-x2)/(f1-f2)+x1;
f2=det(a,b,c,x2)/(x2-rez);
}
cout<<"второе собственное значение: "<<x2<<endl;
cout<<det(a,b,c,x2);
system("pause");
return 0;
}
else
{
xRez1=(-q+sqrt(q*q-4*p*f2))/(2*p);
xRez2=(-q-sqrt(q*q-4*p*f2))/(2*p);

if(abs(xRez2-x2)<abs(xRez1-x2)) xRez1=xRez2;
if(xRez1<x2)x1=xRez1;
else x3=xRez1;
}
}
while(abs(det(a,b,c,xRez1))>eps);
cout<<endl<<setprecision(36)<<"второе собственное значение: "<<xRez1<<endl;
cout<<det(a,b,c,xRez1)<<endl;
system("pause");
return 0;
}
//---------------------------------------------------------------------------


 Оффлайн  Профиль  PM 
 К началу страницы 
+ Ответить 

 Ответить  Открыть новую тему 
1 чел. читают эту тему (гостей: 1, скрытых пользователей: 0)
Пользователей: 0

 





- Текстовая версия 8.09.2025 22:11
500Gb HDD, 6Gb RAM, 2 Cores, 7 EUR в месяц — такие хостинги правда бывают
Связь с администрацией: bu_gen в домене octagram.name