МІНІСТЕРСТВО ОСВІТИ І НАУКИ УКРАЇНИ
СУМСЬКИЙ ДЕРЖАВНИЙ УНІВЕРСИТЕТ
КАФЕДРА ІНФОРМАТИКИ
Курсова робота
з дисципліни «Чисельні методи»
на тему:
«Обчислення характеристичних многочленів, власних значень і власних векторів»
Суми, 2005
Зміст
ЗМІСТ
ТЕОРЕТИЧНІ ДАНІ
ВСТУП
МЕТОД ДАНИЛЕВСЬКА
ВКАЗІВКИ ЩОДО ВИКОРИСТАННЯ ПРОГРАМИ
ПРОГРАМНА РЕАЛІЗАЦІЯ
АНАЛІЗ ПРОГРАМИ
Список використаної літератури
Теоретичні дані
Введення
Велика кількість завдань з механіки, фізики та техніки вимагає знаходження власних значень і власних векторів матриць, тобто таких значень λ, для яких існує нетривіальне рішення однорідної системи лінійних алгебраїчних рівнянь . Тут А-дійсна квадратична матриця порядку n з елементами a jk, а
- Вектор з компонентами x 1, x 2, ..., x n Кожному власному значенню λ i відповідає хоча б одне нетривіальне рішення. Якщо навіть матриця А дійсна, їй власні числа (всі або деякі) і власні вектори можуть бути недійсними. Власні числа є коренями рівняння
, Де Е - одинична матриця порядку n
або
Дане рівняння називається характеристичним рівнянням матриці А. Власним векторах , Яким відповідає власному значенню λ i, називають ненульове рішення однорідної системи рівнянь
. Таким чином, завдання знаходження власних чисел і власних векторів зводиться до знаходження коефіцієнтів характеристичного рівняння, знаходженню його коренів і знаходженню нетривіального рішення системи.
Метод Данилевського
Простий і вишуканий метод знаходження характеристичного многочлена запропонував А. М. Данилевський. Розглянемо ідею методу. Розглянемо матрицю A
Для якої знаходиться характеристичний многочлен, за допомогою подібних перетворень перетвориться до матриці
,
яка має нормальну форму Ф робеніуса, тобто матриця має в явному вигляді в останньому стовпці шукані коефіцієнти характеристичного рівняння. Т. к. подібні матриці мають один і той же характеристичний многочлен, а
, То і
.
Тому для обгрунтування методу досить показати, яким чином з матриці A будується матриця P.
Подібні перетворення матриці A до матриці P відбуваються послідовно. На першому кроці матриця А перетворюється до подібної до неї матриці А (1), в якій передостанній стовпець має необхідний вигляд. На другому кроці матриця А (1) перетворюється на подібну до неї матрицю А (2), в якої вже два передостанніх стовпця мають необхідний вид, і т.д.
На першому кроці матриця А множиться справа на матрицю
і зліва на матрицю їй зворотний
Перший крок дає
,
де
На другому кроці матриця А (1) множиться справа на матрицю
і ліворуч на зворотну до неї матрицю
Очевидно, що елементи матриці
.
Це означає, що два передостанніх стовпця матриці А (2) мають необхідний вид. Продовжуючи цей процес, після n -1 кроків прийдемо до матриці
,
яка має форму Фробеніуса і подібна до вхідних матриці А. При цьому на кожному кроці елементи матриці А (j) знаходяться за елементами матриці А (j -1) також, як ми знаходили елементи матриці А (2) за елементами А (1). При цьому передбачається, що всі елементи відмінні від нуля. Якщо на j-му кроці виявиться, що
, То продовжувати процес у такому вигляді не буде можливо. При цьому можуть виникнути два випадки:
Серед елементів
є хоча б один, відмінний від нуля, наприклад
. Для продовження процесу поміняємо в А (j) місцями перший і
-Го рядка і одночасно 1-й і
, Процес можна продовжувати, тому що стовпці матриці А (j), наведені до необхідного увазі не будуть зіпсовані.
Всі елементи
, Де F - квадратична матриця порядку j, яка має нормальний вигляд Фробеніуса; В-квадратна матриця порядку n - j, але
, Тобто характеристичний многочлен матриці F є дільником характеристичного многочлена матриці А. Для знаходження характеристичного многочлена матриці А необхідно ще знайти характеристичний многочлен матриці В, для якої використовуємо цей же метод.
Підраховано, що кількість операцій множення і ділення, необхідних для отримання характеристичного многочлена матриці порядку n становить n (n -1) (2 n +3) / 2.
На даному етапі роботи ми отримали характеристичний поліном, корінням якого будуть власні числа матриці А. Процедура знаходження коренів полінома n-го ступеня не проста. Тому скористаємося пакетом MathCAD Professional для реалізації даного завдання. Для пошуку коренів звичайного полінома р (х) ступеня n в Mathcad включена дуже зручна функція polyroots (V). Вона повертає вектор всіх коренів многочлена ступеня n, коефіцієнти якого перебувають у векторі V, мають довжину рівну n +1. Зауважимо, що корені полінома можуть бути як речовими, так і комплексними числами. Таким чином ми маємо власні числа, за допомогою яких ми знайдемо власні вектори нашої матриці А. Для знаходження власних векторів скористаємося функцією eigenvec (A, v i), де А-вихідна матриця, v i-власне число, для якого ми шукаємо власний вектор . Ця функція повертає власний вектор дня v i.
Вказівки щодо використання програми
Ця курсова робота виконана на мові програмування 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 ny 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: = ny;
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 ny 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: = ny;
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);
{Процедура додає до i 1 рядок матриці a i 2-у помножену на 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
Власні вектори ,
,
,
Список використаної літератури
Я. М. Григоренко, Н. Д. Панкратова «обчислювальні методи» 1995 р.
В. Д. Гетмнцев «Л інійна алгебра и Лінійне програмування» 2001р.
Д.Мак-Кракена, У. Дорн «Програмування на Фортрані» 1997р.
http://alglib.manual.ru/eigen/danilevsky.php
http://doors.infor.ru/allsrs/alg/index.html
Посилання (links):