МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИУКРАИНЫ
СУМСКИЙГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
КАФЕДРАИНФОРМАТИКИ
Курсовая работа
по дисциплине«Численные методы»
на тему:
«Вычислениехарактеристических многочленов, собственных значений и собственных векторов»
Сумы, 2005
/>/>/>Содержание
СОДЕРЖАНИЕ
ТЕОРЕТИЧЕСКИЕ ДАННЫЕ
ВВЕДЕНИЕ
МЕТОД ДАНИЛЕВСКОГО
УКАЗАНИЯ ПО ПРИМЕНЕНИЮПРОГРАММЫ
ПРОГРАММНАЯ РЕАЛИЗАЦИЯ
АНАЛИЗ ПРОГРАММЫ
СПИСОК ИСПОЛЬЗУЕМОЙЛИТЕРАТУРЫ
/>/>Теоретическиеданные Введение
Большоеколичество задач с механики, физики и техники требует нахождение собственныхзначений и собственных векторов матриц, т.е. таких значений λ, для которыхсуществует нетривиальное решение однородной системы линейных алгебраическихуравнений />.Тут А-действительная квадратичная матрица порядка 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 elsedan(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 in 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 in 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 ij 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
Собственныевекторы />, />,/>,/>
/>/>Список используемой литературы
Я.М.Григоренко,Н.Д.Панкратова «Обчислювальніметоди» 1995р.
В.Д.Гетмнцев «Лінійна алгебра і лінійнепрограмування» 2001р.
Д.Мак-Кракен,У.Дорн «Программирование на ФОРТРАНЕ» 1997г.
alglib.manual.ru/eigen/danilevsky.php
doors.infor.ru/allsrs/alg/index.html