1. ВСТУП
Цілий ряд інженерних задач зводиться до розгляду систем рівнянь, що мають єдине рішення лише в тому випадку, якщо відомо значення деякого вхідного в них параметра. Цей особливий параметр називається характеристичним, або власним, значенням системи. Із завданнями на власні значення інженер стикається в різних ситуаціях. Так, для тензорів напруг власні значення визначають головні нормальні напруження, а власними векторами задаються напрями, пов'язані з цими значеннями. При динамічному аналізі механічних систем власні значення відповідають власним частотам коливань, а власні вектори характеризують моди цих коливань. При розрахунку конструкцій власні значення дозволяють визначати критичні навантаження, перевищення яких призводить до втрати стійкості.
Вибір найбільш ефективного методу визначення власних значень або власних векторів для даної інженерної задачі залежить від ряду факторів, таких, як тип рівнянь, число шуканих власних значень і їх характер. Алгоритми рішення задач на власні значення діляться на дві групи. Ітераційні методи дуже зручні і добре пристосовані для визначення найменшого та найбільшого власних значень. Методи перетворень подібності кілька складніша, зате дозволяють визначити всі власні значення та власні вектори.
У даній роботі будуть розглянуті найбільш поширені методи вирішення задач на власні значення. Однак спочатку наведемо деякі основні відомості з теорії матричного і векторного числень, на яких базуються методи визначення власних значень.
2. ДЕЯКІ ОСНОВНІ ВІДОМОСТІ, НЕОБХІДНІ ПРИ ВИРІШЕННІ ЗАДАЧ НА ВЛАСНІ ЗНАЧЕННЯ
У загальному вигляді завдання на власні значення формулюється таким чином:
A X = l X,
де A - матриця розмірності n х n. Потрібно знайти n скалярних значень l і власні вектори X, відповідні кожному з власних значень.
Основні визначення матричного обчислення
1. Матриця A називається симетричною, якщо
а ij = а ij, де i, j = 1, 2,. . ., N.
Звідси випливає симетрія щодо діагоналі
а kk, де k == 1, 2,. . ., N.
Матриця
1 | 4 | 5 |
4 | 3 | 7 |
5 | 7 | 2 |
є прикладом симетричною.
2. Матриця A називається трехдіагональной, якщо всі її елементи, крім елементів головної і прилеглих до неї діагоналей, дорівнюють нулю. У загальному випадку трехдіагональная матриця має вигляд
* | * | 0 | ||||||
* | * | * | ||||||
* | * | * | ||||||
. | . | . | . | . | . | |||
* | * | * | ||||||
0 | * | * | * | |||||
* | * |
Важливість трехдіагональной форми зумовлена тим, що деякі методи перетворень подібності дозволяють привести довільну матрицю до цього приватного увазі.
3. Матриця A називається ортогональною, якщо
А Т А = Е,
де А т - транспонована матриця A, а Е - одинична матриця. Очевидно, матриця, зворотна ортогональної, еквівалентна транспонованої.
4. Матриці А і В називаються подібними, якщо існує така несінгулярная матриця Р, що справедливе співвідношення
В = Р -1 АР.
Основні властивості власних значень.
1. Всі п власних значень симетричної матриці розмірності ПХП, що складається з дійсних чисел, дійсні. Це корисно пам'ятати, так як матриці, що зустрічаються в інженерних розрахунках, часто бувають симетричними.
2. Якщо власні значення матриці різні, то її власні вектори ортогональні. Сукупність п лінійно незалежних власних векторів утворює базис розглянутого простору. Отже, для сукупності лінійно незалежних власних векторів
X i, де i == 1,. . ., N,
будь-який довільний вектор в тому ж просторі можна виразити через власні вектори. Таким чином,
n
Y = S a i X i.
i = 1
3. Якщо дві матриці подібні, то їх власні значення збігаються. З подоби матриць A і В випливає, що
В = Р -1 АР.
Так як
А Х = l Х,
то
Р -1 АХ = l Р -1 Х.
Якщо взяти Х == Р Y, то
Р -1 АР Y = l Y,
а
У Y == l Y.
Таким чином, матриці A і В не тільки мають однакові власні значення, а й їхні власні вектори зв'язані співвідношенням
Х = Р Y.
4. Помноживши власний вектор матриці на скаляр, одержимо власний вектор тієї ж матриці. Зазвичай всі власні вектори нормують, розділивши кожен елемент власного вектора або на його найбільший елемент, або на суму квадратів всіх інших елементів.
3. Ітераційні методи РІШЕННЯ.
Мабуть, найбільш очевидним способом вирішення задачі на власні значення є їх визначення із системи рівнянь
(A - l E) Х == 0,
яка має ненульове рішення лише у випадку, якщо det (A - l E) = 0. Розкривши визначник, отримаємо многочлен п - й ступеня щодо l, коріння якого і будуть власними значеннями матриці. Для визначення коренів можна скористатися будь-яким з методів, описаних у гл. 2. На жаль, в задачах на власні значення часто зустрічаються кратні корені. Так як ітераційні методи, в цих випадках не гарантують отримання рішення, то для визначення власних значень слід користуватися іншими ітераційними методами.
Визначення найбільшого власного значення методом ітерацій
На рис. 1 показана блок-схема найпростішого ітераційного методу відшукання найбільшого власного значення системи
A Х = l Х.
Процедура починається з пробного нормованого вектора X (0). Цей вектор умножається зліва на матрицю A, і результат прирівнюється твору постійної (власне значення) і нормованого вектору X (0) .. Якщо вектор X (0) збігається з вектором X (0), то рахунок припиняється. В іншому випадку новий нормований вектор використовується як вихідний і вся процедура повторюється. Якщо процес сходиться, то постійний множник відповідає істинному найбільшому власним значенням, а нормований вектор - відповідному власному вектору. Швидкість збіжності цього ітераційного процесу залежить від того наскільки вдало обраний початковий вектор. Якщо він близький до істинного власному вектору, то ітерації сходяться дуже швидко. На швидкість збіжності впливає також і ставлення величин двох найбільших власних значень. Якщо це відношення близьке до одиниці, то збіжність виявляється повільною.
Рис. 1. Блок-схема алгоритму і ітераційного методу розв'язання задач на власні значення.
П Рімера 1
Досліджуємо тривісне напружений стан елемента тіла, представленого на малюнку 2. Матриця напруг для нього має вигляд
10 | 5 | 6 | |
5 | 20 | 4 | * 10 6 Н / м 2 |
6 | 4 | 30 |
Рисунок 2. Т рехосное напружений стан елемента тіла.
Якщо виходити з того, що руйнування відбудеться при максимальній напрузі, то необхідно знати величину найбільшого головного напруги, яка відповідає найбільшому власному значенню матриці напруг. Для знаходження цієї напруги скористаємося методом ітерації Нижче наведена програма для ЕОМ, за допомогою якої ітераційної процедури здійснюється до тих пір, поки різниця між власними значеннями, обчисленими в послідовних ітераціях, не стане менше 0,01%. У програмі використані дві підпрограми - GMPRD з пакету програм для наукових досліджень фірми I ВМ, що служить для перемножування матриць і NORML, нормуються власні вектори за найбільшим елементу.
{************************************************* *********************}
Програма визначення власних значень Програма дозволяє визначить ь найбільше головне напруга (власне значення) для даного тривісного напруженого стану. Застосовується метод ітерацій. Рахунок припиняється, коли зміна власного значення стає менше 0,01 відсотка або число ітерацій перевищує 50.
{************************************************* *********************}
DIMENSION S (3,3), X (3), R (3)
S (1, 1) = 10.E06
S (1,2) = 5.ЕО6
S (2,1) = S (1,2)
S (1,3) = 6.E06
S (3,1) = S (1,3)
S (2,2) = 20.E06
S (2,3) = 4.E06
S (3,2) = S (2,3)
S (3,3) = З0.Е06
X (1) = 1.
Х (2) = 0.0
Х (3) = 0.0
XOLD = 0.0
I = 0
WRITE (6 100)
WRITE (6 101)
WRITE (6 102)
WRITE (6 100)
WRITE (6 104) I, X (1), X (2), X (3)
DO 1 1 = 1,50
CALL GMPRD (S, X, R, 3, 3, 1)
DO 2 J = 1,3
2 X (J) = R (J)
CALL NORML (XLAM, X)
WRITE (6,103) I, XLAM, X (1), X (2), X (3)
IF (ABS ((XOLD-XLAM) / XLAM). LE. 0.0001) GO TO 3
XOLD = XLAM
3 WRITE (6,100)
100 FORMAT (1X 54C'-''))
FORMAT (2X 'ITERATION', зх 'ITERATION', 11X, 'EIGENVECTOR')
FORMAT (3X 'NUMBER ", 6X,' (N / M ** 2) ', 5X,' X (1) ',
6X, 'X (2)', 6X, 'X (3)')
103 FORMAT (1X, I5, 7X, E12.5, 3F10.5)
104 FORMAT (1X, I5, 19X, 3F10.5)
STOP
END
{************************************************* *********************}
SUBROUTINE NORML (XL, X)
DIMENSION X (3)
{************************************************* *********************}
Підпрограма norml.
Ця підпрограма знаходить найбільший з трьох елементів власного вектора і нормує власний вектор з цього найбільшого елементу.
{************************************************* *********************}
# FIND THE LARGEST ELEMENT
XBIG = X (1)
IF (X (2). GT.XBIG) XBIG = X (2)
IF (X (3). GT.XBIG) XBIG = X (3)
# Нормування за XBIG
X (l) = X (1) / XBIG
X (2) = X (2) / XBIG
X (3) = X (3) / XBIG
XL = XBIG
RETURN
END
{************************************************* *********************}
Результат робіт и програми отримуємо у вигляді:
Номер Ітерації | Власне Значення (N / M ** 2) | Власний вектор | ||
X (1) | X (2) | X (3) | ||
0. |
1.00000 | 0. | 0. | ||
0.10000 Е 08 | 1,00000 | 0.50000 | 0.60000 | |
0.26000Е 08 | 0.61923 | 0.66923 | 1.00000 | |
0.36392Е 08 | 0.42697 | 0.56278 | 1.00000 | |
0.34813Е 08 | 0.37583 | 0.49954 | 1.00000 | |
0.34253Е 08 | 0.35781 | 0.46331 | 1.00000 | |
0.34000Е 08 | 0.34984 | 0.44280 | 1.00000 | |
0.33870Е 08 | 0.34580 | 0.43121 | 1.00000 | |
0.33800Е 08 | 0.34362 | 0.42466 | 1.00000 | |
0.33760Е 08 | 0,34240 | 0.42094 | 1.00000 | |
0.33738Е 08 | 0.34171 | 0.41884 | 1.00000 | |
0.33726Е 08 | 0.34132 | 0.41765 | 1.00000 | |
0.33719Е 08 | 0,34110 | 0.41697 | 1.00000 | |
0.33714Е 08 | 0.34093 | 0.41658 | 1.00000 | |
0.33712Е 08 | 0.34091 | 0.41636 | 1.00000 |
Зазначимо, що для досягнення необхідної точності потрібно 14 ітерацій.
Визначення найменшого власного значення методом ітерацій
У деяких випадках доцільно шукати найменше, а не найбільше власне значення. Це можна зробити, попередньо помноживши вихідну систему на матрицю, зворотну A:
А -1 А X = l А -1 X.
Якщо обидві частини цього співвідношення помножимо на 1 / l, то одержимо
1 / l Х = A -1 X.
Ясно, що це вже інше завдання на власне значення, для якої воно дорівнює 1 / l, а розглянутої матрицею є A -1. Максимум 1 / l, досягається при найменшому l. Таким чином, описана вище ітераційної процедури може бути використана для визначення найменшого власного значення нової системи.
Визначення проміжних власних значень методом ітерацій
Знайшовши найбільше власне значення, можна визначити наступне за ним по величині, замінивши вихідну матрицю матрицею, яка містить лише решту власні значення. Використовуємо для цього метод, званий методом вичерпування. Для вихідної симетричної матриці A з відомим найбільшим власним значенням l 1 і власним вектором X 1 можна скористатися принципом ортогональності власних векторів, тобто записати
Х i T Х j = 0 при i <> j і Х i T Х j = 1 при i = j.
Якщо утворити нову матрицю A * відповідно до формули
A * = A-l 1 х 1 Х 1 T,
то її власні значення та власні вектори будуть зв'язані співвідношенням
А * X i = l i X i.
З наведеного вище висловлювання для матриці A * випливає, що
A * Х i = A Х i - l Х 1 Х 1 T X i.
Тут при i = 1 властивість ортогональності дозволяє привести праву частину до виду
A Х 1 - L 1 Х 1.
Але за визначенням власних значень матриці A цей вираз має дорівнювати нулю. Отже, власне значення l 1 Матриці A * дорівнює нулю, а всі інші її власні значення збігаються з власними значеннями матриці A. Таким чином, матриця A * має власні значення 0, l 2, l 3,.. ., L n і відповідні власні вектори Х 1, Х 2, Хз,.. . .... Х n. В результаті виконаних перетворень найбільше власне значення l 1 було вилучено, і тепер, щоб знайти таке найбільше власне значення l 2, можна застосувати до матриці A * звичайний ітераційний метод. Визначивши l 2 і Х 2, повторимо весь процес, використовуючи нову матрицю A **, отриману за допомогою A *, l 2 і Х 2. Хоча на перший погляд здається, що цей процес повинен швидко привести до мети, він має суттєві недоліки. При виконанні кожного кроку похибки у визначенні власних векторів будуть позначатися на точності визначення наступного власного вектора і викликати нагромадження помилок. Тому описаний метод навряд чи застосовний для знаходження більш ніж трьох власних значень, починаючи з найбільшого або найменшого. Якщо потрібно отримати більшу кількість власних значень, слід користуватися методами перетворення подібності.
4. ВИЗНАЧЕННЯ ВЛАСНИХ ЗНАЧЕНЬ МЕТОДАМИ ПЕРЕТВОРЕНЬ ПОДІБНОСТІ
Метод перетворень подібності застосовується з метою отримати з вихідної матриці нову з тими ж власними значеннями, але більш простого вигляду. Очевидно, найкращим спрощенням було б приведення матриці до чисто діагонального вигляду, оскільки в цьому випадку власні значення просто відповідали б елементів матриці, що стоять на головній діагоналі. На жаль, більша частина методів перетворення не дозволяє цього зробити, і доводиться задовольнятися приведенням матриці до трехдіагональной формі.
Метод Якобі
Метод Якобі дозволяє привести матрицю до діагонального вигляду, послідовно, виключаючи всі елементи, що стоять поза головною діагоналі. На жаль, приведення до строго діагонального вигляду вимагає нескінченно великого числа кроків, тому що утворення нового нульового елемента на місці одного з елементів матриці часто веде до появи ненульового елемента там, де раніше був нуль. На практиці метод Якобі розглядають, як итерационную процедуру, яка в принципі дозволяє досить близько підійти до діагональної формі, щоб це перетворення можна було вважати закінченим. У разі симетричної матриці A дійсних чисел перетворення виконується за допомогою ортогональних матриць, отриманих в результаті обертанні в дійсною площині. Обчислення здійснюються таким чином. З вихідної матриці А утворюють матрицю A 1 == Р 1 АР 1 T. При цьому ортогональна матриця Р 1 вибирається так, щоб в матриці А 1 з'явився нульовий елемент, що стоїть поза головної діагоналі. Потім з А 1 за допомогою другої перетворюючої матриці Р 2, утворюють нову матрицю A 2. При цьому Р 2, вибирають так, щоб в A 2 з'явився ще один нульової внедіагональний елемент. Цю процедуру продовжують, прагнучи, щоб на кожному кроці в нуль звертався найбільший внедіагональний елемент. Перетворююча матриця для здійснення зазначеної операції на кожному кроці конструюється таким чином. Якщо елемент а kl матриці А т-1 має максимальну величину, то Рт відповідає
P kk = P ll = cos q,
P kl = - P lk = sin q,
P ii = 1 при i <> k, l, P ij = 0 п ри i <> j.
Матриця А т буде відрізнятися від матриці A m-1 тільки рядками і стовпцями з номерами k і l. Щоб елемент а kl (m) дорівнював нулю, значення q вибирається так, щоб
2 a kl (m-1)
tg 2 q = -------------------------.
a kk (m-1) - a ll (m-1)
k | l | |||||||||||||||||
1 | ||||||||||||||||||