Вычисление характеристических многочленов, собственных значений и собственных векторов
МЕТОД ДАНИЛЕВСКОГО
УКАЗАНИЯ ПО ПРИМЕНЕНИЮ ПРОГРАММЫ
ПРОГРАММНАЯ РЕАЛИЗАЦИЯ
АНАЛИЗ ПРОГРАММЫ
СПИСОК ИСПОЛЬЗУЕМОЙ ЛИТЕРАТУРЫ
Теоретические данные
Введение
Большое количество задач с механики, физики и техники требует нахождение собственных значений и собственных векторов матриц, т.е. таких значений λ, для которых существует нетривиальное решение однородной системы линейных алгебраических уравнений 
. Тут А-действительная квадратичная матрица порядка n с элементами ajk, а 
--вектор с компонентами x1, x2,тАж, xn Каждому собственному значению λi соответствует хотя бы одно нетривиальное решение. Если даже матрица А действительная, ей собственные числа (все или некоторые) и собственные векторы могут быть недействительными. Собственные числа являются корнями уравнения 
, где Е - единичная матрица порядка n
или 
Данное уравнение называется характеристическим уравнением матрицы А. Собственным векторам 
, которым соответствует собственному значению λi, называют ненулевое решение однородной системы уравнений 
. Таким образом, задача нахождения собственных чисел и собственных векторов сводится к нахождению коэффициентов характеристического уравнения, нахождению его корней и нахождению нетривиального решения системы.
Простой и изысканный метод нахождения характеристического многочлена предложил А.М.Данилевский. Рассмотрим идею метода. Рассмотрим матрицу A
Для которой находится характеристический многочлен, при помощи подобных преобразований преобразуется к матрице
,
которая имеет нормальную форму Фробениуса, то есть матрица имеет в явном виде в последнем столбце искомые коэффициенты характеристического уравнения. Т. к. подобные матрицы имеют один и тот же характеристический многочлен, а
, то и 
.
Поэтому для обоснования метода достаточно показать, каким образом из матрицы A строится матрица P.
Подобные преобразования матрицы A к матрице P происходят последовательно. На первом шаге матрица А преобразовывается к подобной до неё матрице А(1), в которой предпоследний столбец имеет необходимый вид. На втором шаге матрица А(1) преобразовывается на подобную к ней матрицу А(2), в которой уже два предпоследних столбца имеют необходимый вид, и т.д.
На первом шаге матрица А умножается справа на матрицу
и слева на матрицу ей обратную
Первый шаг даёт
,
где 
На втором шаге матрица А(1) умножается справа на матрицу
и слева на обратную к ней матрицу
Очевидно, что элементы матрицы
. 
Это означает, что два предпоследних столбца матрицы А(2) имеют необходимый вид. Продолжая этот процесс, после n-1 шагов придем к матрице
,
которая имеет форму Фробениуса и подобная к входной матрице А. При этом на каждом шаге элементы матрицы А(j) находятся по элементам матрицы А(j-1) также, как мы находили элементы матрицы А(2) по элементам А(1). При этом предпологается, что все элементы 
Ваотличные от нуля. Если на j-ом шаге окажется, что 
, то продолжать процесс в таком виде не будет возможно. При этом могут возникнуть два случая:
1. Среди элементов 
Ваесть хотя бы один, отличный от нуля, например 
. Для продолжения процесса поменяем в А(j) местами первый и 
-й строчки и одновременно 1-й и -й столбцы. Такое преобразование матрицы А(j) будет подобным. После того, как получим матрицу 
, процесс можно продолжать, т.к. столбцы матрицы А(j),приведённые к необходимому виду не будут испорчены.
2. Все элементы равны нулю. Тогда матрица А(j) имеет вид 
, где F- квадратичная матрица порядка j, которая имеет нормальный вид Фробениуса; ВтАФквадратная матрица порядка n-j, но 
, то есть характеристический многочлен матрицы F является делителем характеристического многочлена матрицы А. Для нахождения характеристического многочлена матрицы А необходимо еще найти характеристический многочлен матрицы В, для которой используем этот же метод.
Подсчитано, что количество операций умножения и деления, необходимых для получения характеристического многочлена матрицы порядка n составляет n(n-1)(2n+3)/2.
На данном этапе работы мы получили характеристический полином, корнями которого будут собственные числа матрицы А. Процедура нахождения корней полинома n-ой степени не проста. Поэтому воспользуемся пакетом MathCAD Professional для реализации данной задачи. Для поиска корней обычного полинома р(х) степени n в Mathcad включена очень удобная функция polyroots(V). Она возвращает вектор всех корней многочлена степени n, коэффициенты которого находятся в векторе V, имеющим длину равную n+1. Заметим, что корни полинома могут быть как вещественными, так и комплексными числами. Таким образом мы имеем собственные числа, при помощи которых мы найдём собственные векторы нашей матрицы А. Для нахождения собственных векторов воспользуемся функцией eigenvec(A,vi), где А-исходная матрица, vi-собственное число, для которого мы ищем собственный вектор. Данная функция возвращает собственный вектор дня vi.
Указания по применению программы
Данная курсовая работа выполнена на языке программирования Pascal. В курсовую работу входит файл danil.exe. Danil.exe предназначен для нахождения характеристического полинома методом Данилевского. Входными параметрами является размерность матрицы и сама матрица, а выходным тАФ характеристический полином.
Программная реализация
Программный код программы danil.exe
uses wincrt;
label 1;
type mas=array[1.10,1.10]of real;
var A,M,M1,S:mas;
Ваz,max:real;
Ваf,jj,tt,ww,v,h,b,y,i,j,w,k,e,l,q,x,u:byte;
Ваp,o:array[1.10]of real;
Ваt:array [1.10]of boolean;
procedure Umnogenie(b,c:mas; n:byte; var v:mas);
var i,j,k:byte;
begin
for i:=1 to n do
Ваfor j:=1 to n do
Ваbegin
Ваv[i,j]:=0;
Ваfor k:=1 to n do
Ваv[i,j]:=b[i,k]*c[k,j]+v[i,j];
Ваend;
end;
procedure dan(n:byte; var a:mas);
label 1,2;
var y:byte;
begin
For y:=1 to n-1 do
begin
Ваif a[1,n]=0 then
Ваbegin
Ваif y>1 then begin
Ваmax:=abs(a[1,n]);
Ваw:=1;
Ваfor i:=1 to n-y do
Ваif abs(a[i,n])>max then begin max:=abs(a[i,j]); w:=i; end;
Ваif max=0 then
Ваbegin
Ваfor l:=n downto n-y+1 do
Ваbegin
Ваp[f]:=a[l,n];
Ваt[f]:=false;
Ваf:=f-1;
Ваend;
Ваt[f+1]:=true;
Ваx:=x+1;
Ваu:=n-y;
Ваif y=n-1 then begin o[q]:=a[1,1]; q:=q+1; end else dan(u,a);
Ваgoto 2;
Ваend;
Ваfor j:=1 to n do
Ваbegin
Ваz:=a[1,j];
Ваa[1,j]:=a[w,j];
Ваa[w,j]:=z;
Ваend;
Ваfor k:=1 to n do
Ваbegin
Ваz:=a[k,1];
Ваa[k,1]:=a[k,w];
Ваa[k,w]:=z;
Ваend;
Ваgoto 1;
Ваend
Ваelse
Ваbegin
Ваmax:=abs(a[1,2]);
Ваw:=1;e:=2;
Ваfor i:=1 to n-1 do
Ваif abs(a[i,n])>max then begin max:=abs(a[i,j]); w:=i; e:=n; end;
Ваfor j:=2 to n do
Ваif abs(a[1,j])>max then begin max:=abs(a[i,j]); w:=1; e:=j; end;
Ваif abs(a[n,1])>max then begin max:=abs(a[n,1]); w:=n; e:=1; end;
Ваif max=0 then
Ваbegin
Ваo[q]:=a[n,n];
Ваq:=q+1;
Ваu:=n-1;
Ваif n=2 then begin o[q]:=a[1,1]; q:=q+1; o[q]:=a[n,n]; q:=q+1; end else dan(u,a);
Ваgoto 2;
Ваend;
Ваif (w>1) and (e=n) then
Ваbegin
Ваfor j:=1 to n do
Ваbegin
Ваz:=a[1,j];
Ваa[1,j]:=a[w,j];
Ваa[w,j]:=z;
Ваend;
Ваfor k:=1 to n do
Ваbegin
Ваz:=a[k,1];
Ваa[k,1]:=a[k,w];
Ваa[k,w]:=z;
Ваend;
Ваgoto 1;
Ваend;
Ваif (w=n) and (e=1) then
Ваbegin
Ваfor j:=1 to n do
Ваbegin
Ваz:=a[1,j];
Ваa[1,j]:=a[n,j];
Ваa[n,j]:=z;
Ваend;
Ваfor k:=1 to n do
Ваbegin
Ваz:=a[k,1];
Ваa[k,1]:=a[k,n];
Ваa[k,n]:=z;
Ваend;
Ваgoto 1;
Ваend;
Ваif w=1 then
Ваbegin
Ваfor j:=1 to n do
Ваbegin
Ваz:=a[n,j];
Ваa[n,j]:=a[e,j];
Ваa[e,j]:=z;
Ваend;
Ваfor k:=1 to n do
Ваbegin
Ваz:=a[k,n];
Ваa[k,n]:=a[k,e];
Ваa[k,e]:=z;
Ваend;
Ваgoto 1;
Ваend;
Ваend;
end;
1:
Ваfor i:=1 to n do
Ваfor j:=1 to n do
Ваif i<>(j+1) then M[i,j]:=0
Ваelse M[i,j]:=1;
Ваfor i:=1 to n do
Ваfor j:=1 to n do
Ваif (i+1)<>j then M1[i,j]:=0
Ваelse M1[i,j]:=1;
Ваfor i:=1 to n do
Ваif i<>n then begin M[i,n]:=a[i,n]; M1[i,1]:=-a[i+1,n]/a[1,n]; end
Ваelse begin M[i,n]:=a[i,n]; M1[i,1]:=1/a[1,n]; end;
ВаUmnogenie(M1,A,n,S);
ВаUmnogenie(S,M,n,A);
if y=n-1 then
begin
Ваfor l:=n downto 1 do
Ваbegin
Ваp[f]:=a[l,n];
Ваt[f]:=false;
Ваf:=f-1;
Ваend;
Ваt[f+1]:=true;
Ваx:=x+1;
end;
end;
2:
end;
begin
writeln('Vvedite razmernost` matrici A');
readln(ww);
f:=ww;
for i:=1 to ww do
begin
Ваfor j:=1 to ww do
Ваbegin
Ваwrite('a[',i,j,']=');
ВаReadln(A[i,j]);
Ваend;
end;
q:=1;
x:=0;
dan(ww,a);
for i:=1 to q-1 do
writeln('Koren` har-ogo ur-iya=',o[i]:2:2);
writeln;
i:=ww+1;
if (x=1)or(x>1) then
Ваbegin
Ваfor v:=1 to x do
Ваbegin
Ваtt:=0;
Ваrepeat
Ваtt:=tt+1;
Ваi:=i-1;
Ваuntil t[i]<>false;
Ваwrite('l^',tt,' + ');
Ваfor jj:=ww downto i do
Ваbegin
Ваtt:=tt-1;
Ваwrite(-p[jj]:2:2,'*l^',tt,' + ');
Ваend;
Ваww:=i-1;
Ваwriteln;
Ваend;
Ваend;
end.
Получение формы Жордано: form.exe
uses wincrt;
label 1;
type mas=array[1.10,1.10]of real;
var A,M,M1,S,R,R1,A1:mas;
Ваz,max:real;
Ваf,jj,tt,ww,v,h,b,y,i,j,w,k,e,l,q,x,u,n1:byte;
Ваp,o:array[1.10]of real;
Ваt:array [1.10]of boolean;
procedure Umnogenie(b,c:mas; n:byte; var v:mas);
var i,j,k:byte;
begin
for i:=1 to n do
Ваfor j:=1 to n do
Ваbegin
Ваv[i,j]:=0;
Ваfor k:=1 to n do
Ваv[i,j]:=b[i,k]*c[k,j]+v[i,j];
Ваend;
end;
procedure dan(n:byte; var a:mas);
label 1,2;
var y:byte;
begin
For y:=1 to n-1 do
begin
Ваif a[1,n]=0 then
Ваbegin
Ваif y>1 then begin
Ваmax:=abs(a[1,n]);
Ваw:=1;
Ваfor i:=1 to n-y do
Ваif abs(a[i,n])>max then begin max:=abs(a[i,j]); w:=i; end;
Ваif max=0 then
Ваbegin
Ваfor l:=n downto n-y+1 do
Ваbegin
Ваp[f]:=a[l,n];
Ваt[f]:=false;
Ваf:=f-1;
Ваend;
Ваt[f+1]:=true;
Ваx:=x+1;
Ваu:=n-y;
Ваif y=n-1 then begin o[q]:=a[1,1]; q:=q+1; end else dan(u,a);
Ваgoto 2;
Ваend;
Ваfor j:=1 to n do
Ваbegin
Ваz:=a[1,j];
Ваa[1,j]:=a[w,j];
Ваa[w,j]:=z;
Ваend;
Ваfor k:=1 to n do
Ваbegin
Ваz:=a[k,1];
Ваa[k,1]:=a[k,w];
Ваa[k,w]:=z;
Ваend;
Ваgoto 1;
Ваend
Ваelse
Ваbegin
Ваmax:=abs(a[1,2]);
Ваw:=1;e:=2;
Ваfor i:=1 to n-1 do
Ваif abs(a[i,n])>max then begin max:=abs(a[i,j]); w:=i; e:=n; end;
Ваfor j:=2 to n do
Ваif abs(a[1,j])>max then begin max:=abs(a[i,j]); w:=1; e:=j; end;
Ваif abs(a[n,1])>max then begin max:=abs(a[n,1]); w:=n; e:=1; end;
Ваif max=0 then
Ваbegin
Ваo[q]:=a[n,n];
Ваq:=q+1;
Ваu:=n-1;
Ваif n=2 then begin o[q]:=a[1,1]; q:=q+1; o[q]:=a[n,n]; q:=q+1; end else dan(u,a);
Ваgoto 2;
Ваend;
Ваif (w>1) and (e=n) then
Ваbegin
Ваfor j:=1 to n do
Ваbegin
Ваz:=a[1,j];
Ваa[1,j]:=a[w,j];
Ваa[w,j]:=z;
Ваend;
Ваfor k:=1 to n do
Ваbegin
Ваz:=a[k,1];
Ваa[k,1]:=a[k,w];
Ваa[k,w]:=z;
Ваend;
Ваgoto 1;
Ваend;
Ваif (w=n) and (e=1) then
Ваbegin
Ваfor j:=1 to n do
Ваbegin
Ваz:=a[1,j];
Ваa[1,j]:=a[n,j];
Ваa[n,j]:=z;
Ваend;
Ваfor k:=1 to n do
Ваbegin
Ваz:=a[k,1];
Ваa[k,1]:=a[k,n];
Ваa[k,n]:=z;
Ваend;
Ваgoto 1;
Ваend;
Ваif w=1 then
Ваbegin
Ваfor j:=1 to n do
Ваbegin
Ваz:=a[n,j];
Ваa[n,j]:=a[e,j];
Ваa[e,j]:=z;
Ваend;
Ваfor k:=1 to n do
Ваbegin
Ваz:=a[k,n];
Ваa[k,n]:=a[k,e];
Ваa[k,e]:=z;
Ваend;
Ваgoto 1;
Ваend;
Ваend;
end;
1:
Ваfor i:=1 to n do
Ваfor j:=1 to n do
Ваif i<>(j+1) then M[i,j]:=0
Ваelse M[i,j]:=1;
Ваfor i:=1 to n do
Ваfor j:=1 to n do
Ваif (i+1)<>j then M1[i,j]:=0
Ваelse M1[i,j]:=1;
Ваfor i:=1 to n do
Ваif i<>n then begin M[i,n]:=a[i,n]; M1[i,1]:=-a[i+1,n]/a[1,n]; end
Ваelse begin M[i,n]:=a[i,n]; M1[i,1]:=1/a[1,n]; end;
ВаUmnogenie(M1,A,n,S);
ВаUmnogenie(S,M,n,A);
if y=n-1 then
begin
Ваfor l:=n downto 1 do
Ваbegin
Ваp[f]:=a[l,n];
Ваt[f]:=false;
Ваf:=f-1;
Ваend;
Ваt[f+1]:=true;
Ваx:=x+1;
end;
end;
2:
end;
procedure ObrMatr(A:mas;Var AO:mas; n:byte);
Ваconst e=0.00001;
Ваvar i,j:integer;
Ваa0:mas;
Ваprocedure MultString(var A,AO:mas;i1:integer;r:real);
Ваvar j:integer;
Ваbegin
Ваfor j:=1 to n do
Ваbegin
ВаA[i1,j]:=A[i1,j]*r;
ВаAO[i1,j]:=AO[i1,j]*r;
Ваend;
Ваend;
Ваprocedure AddStrings(var A,AO:mas;i1,i2:integer;r:real);
Ва{Процедура прибавляет к i1 строке матрицы a i2-ю умноженную на r}
Ваvar j:integer;
Ваbegin
Ваfor j:=1 to n do
Ваbegin
ВаA[i1,j]:=A[i1,j]+r*A[i2,j];
ВаAO[i1,j]:=AO[i1,j]+r*AO[i2,j];
Ваend;
Ваend;
Ваfunction Sign(r:real):shortint;
Ваbegin
Ваif (r>=0) then sign:=1
Ваelse sign:=-1;
Ваend;
Ваbegin {начало основной процедуры}
Ваfor i:=1 to n do
Ваfor j:=1 to n do
Ваa0[i,j]:=A[i,j];
Ваfor i:=1 to n do
Ваbegin {К i-той строке прибавляем (или вычитаем)
Ваj-тую строку взятую со знаком i-того
Ваэлемента j-той строки. Таким образом,
Вана месте элемента a[i,i] возникает сумма
Вамодулей элементов i-того столбца (ниже i-той строки)
Вавзятая со знаком бывшего элемента a[i,i],
Варавенство нулю которой говорит о несуществовании
Ваобратной матрицы }
Ваfor j:=i+1 to n do
ВаAddStrings(A,AO,i,j,sign(A[i,i])*sign(A[j,i])); { Прямой ход }
Ваif (abs(A[i,i])>e) then
Ваbegin
ВаMultString(a,AO,i,1/A[i,i]);
Ваfor j:=i+1 to n do
ВаAddStrings(a,AO,j,i,-A[j,i]);
Ваend
Ваelse begin writeln('Обратной матрицы не существует.');
Ваhalt;
Ваend
Ваend;{Обратный ход:}
Ваif (A[n,n]>e) then begin
Ваfor i:=n downto 1 do
Ваfor j:=1 to i-1 do
Ваbegin
ВаAddStrings(A,AO,j,i,-A[j,i]);
Ваend; end
Ваelse writeln('Обратной матрицы не существует.');
Ваend;
procedure EdMatr(Var E:mas; n:byte);
Ваvar i,j:byte;
Ваbegin
Ваfor i:=1 to n do
Ваfor j:=1 to n do
Ваif i<>j then E[i,j]:=0 else E[i,i]:=1;
Ваend;
{procedure UmnogMatr(A,F:mas; Var R:mas; n:byte);
ВаVar s:real;
Ваl,i,j:byte;
Ваbegin
Ваfor i:=1 to n do
Ваfor j:=1 to n do
Ваbegin
Ваs:=0;
Ваfor l:=1 to n do
Ваs:=s+A[i,l]*F[l,j];
ВаR[i,j]:=s;
Ваend;
Ваend; }
begin
writeln('Vvedite razmernost` matrici A');
readln(ww);
f:=ww;
n1:=ww;
for i:=1 to ww do
begin
Ваfor j:=1 to ww do
Ваbegin
Ваwrite('a[',i,j,']=');
ВаReadln(A[i,j]);
ВаA1[i,j]:=A[i,j];
Ваend;
end;
q:=1;
x:=0;
dan(ww,a);
for i:=1 to q-1 do
writeln('Koren` har-ogo ur-iya=',o[i]:2:2);
writeln;
i:=ww+1;
if (x=1)or(x>1) then
Ваbegin
Ваfor v:=1 to x do
Ваbegin
Ваtt:=0;
Ваrepeat
Ваtt:=tt+1;
Ваi:=i-1;
Ваuntil t[i]<>false;
Ваwrite('l^',tt,' + ');
Ваfor jj:=ww downto i do
Ваbegin
Ваtt:=tt-1;
Ваwrite(-p[jj]:2:2,'*l^',tt,' + ');
Ваend;
Ваww:=i-1;
Ваwriteln;
Ваend;
Ваend;
for i:=1 to n1 do
Ваbegin
Ваfor j:=1 to n1 do
Ваread(R[i,j]);
Ваreadln;
Ваend;
EdMatr(R1,n1);
ObrMatr(R,R1,n1);
Umnogenie(R1,A1,n1,A);
Umnogenie(A,R,n1,M1);
for i:=1 to n1 do
Ваbegin
Ваfor j:=1 to n1 do
Ваwrite(' ',M1[i,j]:2:3,' ');
Ваwriteln;
Ваend;
end.
Протестируем работу программы на примере. Пусть имеем матрицу А
Характеристический полином имеет вид:
Собственные числа 20.713, 4.545, 2.556, -5.814
Собственные векторы 
, 
,
,
Список используемой литературыЯ.М.Григоренко, Н.Д.Панкратова ВлОбчислювальнi методиВ» 1995р.
В.Д.Гетмнцев ВлЛiнiйна алгебра i лiнiйне програмуванняВ» 2001р.
Д.Мак-Кракен, У.Дорн ВлПрограммирование на ФОРТРАНЕВ» 1997г.
http://alglib.manual.ru/eigen/danilevsky.php
http://doors.infor.ru/allsrs/alg/index.html
Вместе с этим смотрят:
РЖнварiантнi пiдпростори. Власнi вектори i власнi значення лiнiйного оператора
РЖнтегральнi характеристики векторних полiв
РЖнтерполювання функцiй
Автокорреляционная функция. Примеры расчётов
Аксонометричнi проекцii