Pobierz Wartości i wektory własne - charakterystyka i więcej Skrypty w PDF z Engineering tylko na Docsity!
WydziaPolitechnika Krakowskał Inżynierii Lądowej
Wartości i wektory własne
Wektorem własnym macierzy A [ (^) n × n ]nazywamy każdy niezerowy wektor V , który zachowuje kie- runek po wykonaniu mnożenia przez tę macierz: A ⋅ V = λ ⋅ V. (1) Wielkość λ jest wartością własną macierzy A odpowiadającą wektorowi własnemu V. Rzeczywista i symetryczna macierz A posiada wyłącznie rzeczywiste wartości własne, oraz zbiór n liniowo niezależnych wektorów własnych oznaczonych V 1 , V 2 ,..., Vn. Wektory te oczywiście nie są wyznaczone jednoznacznie, ponieważ każdy z nich może zostać pomnożony przez niezerową stałą, dalej pozostając wektorem własnym (definicja).
Każdy wektor własny ma odpowiadającą wartość własną, oznaczoną λ 1 , λ 2 ,..., λ n. Wartości własne
dla danej macierzy są określone jednoznacznie. Zbiór wszystkich wartości własnych macierzy na- zywamy spektrum tej macierzy. Równaniem charakterystycznym macierzy A nazywamy równanie:
det ( A − λ ⋅ I ) = p ( λ) = 0 , (2)
w którym I (^) [ n × n ] oznacza macierz jednostkową. Wyrażenie stojące po lewej stronie tego równania jest wielomianem p stopnia n zmiennej λ. Tak więc wyznaczenie spektrum macierzy A jest równoznaczne z wyznaczeniem wszystkich pierwiastków równania (2). Macierz A jest pierwiastkiem własnego równania charakterystycznego, czyli:
p ( A ) = 0. (3)
Twierdzenie powyższe nosi nazwę twierdzenia Cayleya – Hamiltona. Ponadto, jeżeli q ( ) x jest wie- lomianem a λ wartością własną macierzy A to q (λ )jest wartością własną macierzy q ( A ).
Przekształcenie przez podobieństwo nie zmienia wartości własnych macierzy, czyli jeżeli istnieje odwracalna macierz P [ (^) n × n ]to wartości własne macierzy P −^1 ⋅ A ⋅ P i macierzy A są równe.
Przekształcenie ortogonalne nie zmienia wartości własnych macierzy, czyli jeżeli istnieje odwracal- na macierz Q [ (^) n × n ]taka, że Q T^ ⋅ Q = I to wartości własne macierzy Q T^ ⋅ A ⋅ Q i macierzy A są równe (wniosek z twierdzenia poprzedniego). Dla symetrycznej dodatnio określonej macierzy A i niezerowego wektora X zachodzi związek:
0 λmin ⋅ ≤ λmax
X X
X A X
T
T
. (4)
Do oszacowania spektrum wartości własnych dowolnej macierzy kwadratowej możemy użyć twier- dzenia Gerszgorina, natomiast w przypadku symetrycznej macierzy trójdiagonalnej możemy zasto- sować ciągi Sturma do określenia z dowolną dokładnością przedziałów, w których znajdują się jej poszczególne wartości własne. Numeryczne wyznaczanie wartości własnych macierzy jest jednym z bardziej żmudnych zadań podstawowej analizy numerycznej. Jednak ze względu na zastosowania praktyczne zostało opraco- wanych wiele metod służących do rozwiązania tego zagadnienia. Metody te można podzielić na trzy kategorie:
WydziaPolitechnika Krakowskał Inżynierii Lądowej
- metody pozwalające na znalezienie wszystkich wartości i wektorów własnych, np. metody: Jacobiego zwana również metodą transformacji ortogonalnych, Lanczosa i Hausholdera w połączenu z metodą rozkładu Q ⋅ R ;
- metody pozwalające na wyznaczenie wybranej wartości i odpowiadającego jej wektora wła- snego, np. metoda potęgowa, w ogólności pozwalająca na wyznaczenie wartości własnej najbliższej podanej liczbie i odpowiadającego jej wektora własnego;
- metody pozwalające na obliczenie grupy wartości własnych i odpowiadających im wekto- rów własnych. Przed wyborem metody obliczeń należy wobec powyższego, ze względu na efektywność czasową obliczeń, zastanowić się, czy będą nam potrzebne wszystkie wartości własne analizowanej macie- rzy, czy tylko ich wybrany podzbiór, czy też jedna jedyna wartość własna spełniająca określone wa- runki. W dalszej części pracy zajmiemy się dwiema metodami wyznaczania wartości własnych rzeczywi- stej macierzy symetrycznej. Metoda Jacobiego Metoda Jacobiego (transformacji ortogonalnych) polega na wykonaniu na wyjściowej macierzy A [ (^) n × n ] ciągu transformacji ortogonalnych, w wyniku których macierz ta zostanie doprowadzona do postaci diagonalnej D [ (^) n × n ]. W macierzy diagonalnej na przekątnej głównej znajdą się wartości wła- sne macierzy wyjściowej, natomiast wektory własne odpowiadające tym wartościom własnym będą zapisane w kolumnach macierzy W [ (^) n × n ].
A = W ⋅ D ⋅ W^ T , (5) gdzie:
W = Q^ { m^ }^^ ⋅... ⋅ Q { } i^ ⋅...⋅ Q { }^2 ⋅ Q { }^1. (6)
jest iloczynem macierzy Q { } i definiujących kolejne transformacje ortogonalne i = 1 , 2 ,..., m.
Transformacje ortogonalne są wykonywane w taki sposób, aby kolejno zerować pary elementów macierzy A rozmieszczonych symetrycznie względem przekątnej głównej i największych co do modułu. Tak więc algorytm metody można rozbić na następujące kroki:
- Wybór elementu wiodącego transformacji ortogonalnej, czyli wyznaczenie indeksów p i q elementu macierzy A największego co do modułu, a nie leżącego na przekątnej głównej tej ma- cierzy: { } { } kli klkl n pq apqi a = ≠ , : = (^) , max 1 , 2 ,...,. (7)
- Wyznaczenie macierzy transformacji Q {^ i }tak, aby elementy nowej macierzy:
A^ {^ i^ +^1 }^ = Q { } iT^ ⋅ A { } i^ ⋅ Q { } i (8) spełniały warunek: a^ { pq^ i +^1 }^ = a { qp^ i +^1 }^ = 0 , (9) i wyznaczenie elementów nowej macierzy A { i^ +^1 } .Jak łatwo sprawdzić, jeżeli macierz Q {^ i } przyjmiemy jako macierz jednostkową, w której zmieniono jedynie elementy:
WydziaPolitechnika Krakowskał Inżynierii Lądowej
t^2 + 2 ⋅ η ⋅ t − 1 = 0 , (18)
czyli zwykłe równanie kwadratowe ze względu na pomocniczą zmienną t. Pierwiastki tego równania można wyznaczyć posługując się dobrze znanymi wzorami, jednak, ze względu na po- prawę stabilności numerycznej rozwiązania, do wyznaczenia mniejszego z nich korzystniej jest stosować formułę:
sgn
t η. (19)
Gdyby jednakże η było tak duże, że η 2 powodowałoby wystąpienie błędu nadmiaru w trakcie
wykonywania obliczeń, należy zastosować formułę:
t^1. (20)
Po wyznaczeniu wartości zmiennej pomocniczej t powracamy do wyjściowych niewiadomych c i s , korzystając z zależności (12) i (17):
s t c
t
c = ⋅
- Obliczenie macierzy wektorów własnych na podstawie zależności:
W^ { i^^ +^1 }^ = W { } i^ ⋅ Q { } i , (22) przy czym początkową macierzą W {^0 } jest macierz jednostkowa. Stosowne wzory transforma- cyjne przyjmą postać: { } { } { } w {^ }^ s w { }^ c w { }^ r n
w c w s w rqi rpi rqi
rpi rpi rqi 1 1 ,^2 ,...,
1
, (23)
która wynika z faktu, że w kolejnych iteracjach zmieniają się wyłącznie kolumny p i q macie- rzy W.
- Sprawdzenie warunku zakończenia obliczeń: { }
1
1 i
i sn
sp (^) , (24)
w którym licznik i mianownik ułamka są odpowiednio równe: { } { }
{ } { 1 } 1 , 2 ,...,
1
1 , 1 , 2 ,...,
1
max
max
=
=≠
=
ii^ i i n
i
iji iji j n
i
sn a
sp a (25)
dla normy maksimum, lub:
=
=≠ =
−
n i
i n iii
i j
n i
n j iji n n
i
sn a
sp a
1
1 1 1 2
1 1
1 1 1 2 2 (26)
WydziaPolitechnika Krakowskał Inżynierii Lądowej
dla normy średniokwadratowej. Jeżeli warunek (24) jest spełniony dla narzuconej przez użytkownika dokładności ε , to kończymy obliczenia, w przeciwnym przypadku kontynuujemy pracę poczynając od wyznaczenia nowego elementu wiodącego a { (^) pqi^ +^1 }(7). Tak obrana metoda postępowania będzie zbieżna, ponieważ jeśli ob- liczymy sumę kwadratów elementów macierzy A { i^ }leżących poza przekątną główną i oznaczymy ją przez S { } i :
i j
n i
n j
Si^ aiji = (^) ≠=
to wprost z (14) wynika, że:
S^ { i^^ +^1 }^ = S { } i^ − 2 ⋅( a { } pqi )^2 , (28)
czyli, że licznik ułamka we wzorze (24) będzie monotonicznie malał do zera w kolejnych itera- cjach, a więc dla dowolnie małego dodatniego ε warunek (24) zostanie spełniony po skończonej liczbie iteracji. Metoda Jacobiego jest efektywna dla macierzy A o niewielkich rozmiarach (rzędu 10). Do prowa- dzenia obliczeń na dużych macierzach należy stosować inne, bardziej efektywne metody, np. meto- dę Hausholdera lub Lanczosa aby doprowadzić macierz wyjściową do postaci trójdiagonalnej a na- stępnie metodę rozkładu Q ⋅ R do wyznaczenia wartości i wektorów własnych tak zredukowanej macierzy. Algorytm postępowania przedstawia się następująco:
1. Wybrać element wiodący (czyli wyznaczyć p i q ) (7);
2. Obliczyć η i t (16), (19);
3. Obliczyć c i s (21); 4. Wyznaczyć A { i^ +^1 }(8) lub (13), (14); 5. Wyznaczyć W { i^ +^1 }(22) lub (23); 6. Sprawdzić warunek (24). W zależności od wyniku sprawdzenia albo należy wrócić do 1., albo zakończyć pracę. Dla lepszej ilustracji sposobu postępowania przedstawimy go na przykładzie obliczania wszystkich wartości własnych i odpowiadających im wektorów własnych macierzy A [ (^) 4 × 4 ] z dokładnością ε = 0 , 0001 :
{ } ⎥
− −
−
− − −
−
4 1 5 7
1 2 2 5
2 6 2 1
3 2 1 4 A A^0. (29)
WydziaPolitechnika Krakowskał Inżynierii Lądowej
Iteracja druga:
- Element wiodący (największy co do modułu element macierzy A {^1 } spoza przekątnej głównej) p = 1 , q = 3.
- Współczynniki η i t : { } { } { } { } { } { }
( )
sgn( 0 , 323308 ) 1
sgn( )
2 12 ,^0901702 , 9535753 0 ,^323308
2 2
131
331 111
111 131 331
η η
t η
a η a a
a a a
. (37)
- Współczynniki c i s :
( ) 0 , 727657 0 , 808588 0 , 588375
2 2 = ⋅ =− ⋅ = −
s t c
t
c (^). (38)
- Macierz A { }^2 :
{ } ⎥
−
0 0 , 000 0 , 000 1 , 000
0 , 588 0 , 000 0 , 809 0 , 000
0 , 000 1 , 000 0 , 000 0 , 000
0 , 809 0 , 000 0 , 588 0 , 000 1
0 0 0 1
s c
c s Q. (39)
A^ { }^2 = Q { }^1 T^ ⋅ A { }^1 ⋅ Q { }^1 =
− − −
− −
− − −
− − − −
− − −
− −
2 , 326 1 , 902 1 , 693 10 , 090
0 , 000 2 , 127 1 , 059 1 , 693
0 , 926 6 , 000 2 , 127 1 , 902
5 , 149 0 , 926 0 , 000 2 , 326
0 , 0000 , 0000 , 0001 , 000
0 , 5880 , 0000 , 8090 , 000
0 , 0001 , 0000 , 0000 , 000
0 , 8090 , 0000 , 5880 , 000
2 , 877 1 , 9020 , 000 10 , 090
2 , 954 1 , 1761 , 090 0 , 000
2 , 000 6 , 0001 , 176 1 , 902
3 , 000 2 , 0002 , 954 2 , 877
0 , 0000 , 000 0 , 0001 , 000
0 , 5880 , 000 0 , 8090 , 000
0 , 0001 , 000 0 , 0000 , 000
0 , 8090 , 000 0 , 5880 , 000
. (40)
- Macierz W { }^2 :
W^ {^2 }^ = W { }^1 ⋅ Q { }^1 =
− − − 0 , 309 0 , 000 0 , 425 0 , 851
0 , 501 0 , 000 0 , 688 0 , 526
0 , 000 1 , 000 0 , 000 0 , 000
0 , 809 0 , 000 0 , 588 0 , 000
0 , 000 0 , 000 0 , 000 1 , 000
0 , 588 0 , 000 0 , 809 0 , 000
0 , 000 1 , 000 0 , 000 0 , 000
0 , 809 0 , 000 0 , 588 0 , 000
0 , 000 0 , 000 0 , 526 0 , 851
0 , 000 0 , 000 0 , 851 0 , 526
0 , 000 1 , 000 0 , 000 0 , 000
1 , 000 0 , 000 0 , 000 0 , 000
. (41)
- Warunek zakończenia obliczeń (w normie maksimum): { } { }
{ } (^) max { } (^10) , 090170
max 2 , 326205
2 1 , 2 , 3 , 4
2
2 , 1 , 2 , 3 , 4
2
=
=≠
i^ ii
i j ij ij
sn a
sp a , (42)
{ } sn { }^ ε
sp (^) = = 0 , 230542 > 10 , 090170
2
2 , (43)
czyli dalej konieczne jest wykonanie następnej iteracji.
WydziaPolitechnika Krakowskał Inżynierii Lądowej
Iteracja trzecia:
- Element wiodący (największy co do modułu element macierzy A {^2 } spoza przekątnej głównej) p = 1 , q = 4.
- Współczynniki η i t : { } { } { } { } { } { }
( )
sgn( 3 , 275584 ) 1
sgn( )
2 10 ,^09017022 , 3262055 ,^1491903 ,^275584
2 2
142
442 112
112 142 442
η η
t η
a η a a
a a a
. (44)
- Współczynniki c i s :
( ) 0 , 1492450 , 989046 0 , 147610
2 2 = ⋅ =− ⋅ = −
s t c
t
c (^). (45)
- Macierz A { }^3 :
{ } ⎥
−
0 , 148 0 , 000 0 , 000 0 , 989
0 , 000 0 , 000 0 , 000 0 , 000
0 , 000 1 , 000 0 , 000 0 , 000
0 , 989 0 , 000 0 , 000 0 , 148 2
0 0
s c
c s Q. (46)
A^ { }^3 = Q { }^2 T^ ⋅ A { }^2 ⋅ Q { }^2 =
− − −
− − −
− − −
− −
− − − −
− −
− − −
− −
0 , 000 1 , 745 1 , 674 10 , 437
0 , 250 2 , 127 1 , 059 1 , 674
1 , 196 6 , 000 2 , 127 1 , 745
5 , 496 1 , 196 0 , 250 0 , 000
0 , 1480 , 0000 , 0000 , 989
0 , 0000 , 0001 , 0000 , 000
0 , 0001 , 0000 , 0000 , 000
0 , 9890 , 0000 , 0000 , 148
2 , 326 1 , 902 1 , 693 10 , 090
0 , 000 2 , 127 1 , 059 1 , 693
0 , 926 6 , 000 2 , 127 1 , 902
5 , 149 0 , 926 0 , 000 2 , 326
0 , 1480 , 0000 , 000 0 , 989
0 , 0000 , 0001 , 000 0 , 000
0 , 0001 , 0000 , 000 0 , 000
0 , 9890 , 0000 , 000 0 , 148
. (47)
- Macierz W { }^3 :
W^ {^3 }^ = W { }^2 ⋅ Q { }^2 =
− − − −
−
0 , 432 0 , 000 0 , 425 0 , 796
0 , 417 0 , 000 0 , 688 0 , 594
0 , 000 1 , 000 0 , 000 0 , 000
0 , 800 0 , 000 0 , 588 0 , 119
0 , 148 0 , 000 0 , 000 0 , 989
0 , 000 0 , 000 1 , 000 0 , 000
0 , 000 1 , 000 0 , 000 0 , 000
0 , 989 0 , 000 0 , 000 0 , 148
0 , 309 0 , 000 0 , 425 0 , 851
0 , 501 0 , 000 0 , 688 0 , 526
0 , 000 1 , 000 0 , 000 0 , 000
0 , 809 0 , 000 0 , 588 0 , 000
. (48)
- Warunek zakończenia obliczeń (w normie maksimum): { } { }
{ } (^) max { } (^10) , 437343
max 2 , 127302
3 1 , 2 , 3 , 4
3
3 , 1 , 2 , 3 , 4
3
=
=≠
i^ ii
i j ij ij
sn a
sp a , (49)
{ } sn { }^ ε
sp (^) = = 0 , 203816 > 10 , 437343
3
3 , (50)
czyli dalej konieczne jest wykonanie następnej iteracji.
WydziaPolitechnika Krakowskał Inżynierii Lądowej
Metoda potęgowa Metoda potęgowa w wersji podstawowej służy do wyznaczenia maksymalnej co do modułu warto- ści własnej i odpowiadającego jej wektora własnego. Jeżeli jednak wykorzystamy fakty, że wartości własne macierzy odwrotnej są odwrotnościami wartości własnych macierzy danej oraz że modyfi- kacja macierzy polegająca na dodaniu do elementów jej przekątnej głównej pewnej liczby d po- woduje przesunięcie wszystkich wartości własnych tej macierzy o tę samą liczbę d ( przesunięcie widma ), możemy przy pomocy metody potęgowej wyznaczyć wartość własną najbliższą dowolnej liczbie. Metoda potęgowa polega na wykonaniu ciągu mnożeń przyjętego wektora startowego przez ma- cierz, której dominującej wartości własnej i odpowiadającego wektora własnego poszukujemy. Pro- cedura taka jest równoznaczna pomnożeniu wektora startowego przez macierz wyjściową podnie- sioną do potęgi równej liczbie iteracji. Wektor będący iloczynem zmierza wtedy do wektora odpo- wiadającego największej wartości własnej, samą wartość własną możemy natomiast obliczyć z wy- rażenia zwanego ilorazem Rayleigha (4). Uzasadnienie poprawności metody jest następujące. Roz- ważmy dowolny wektor X { }^0 oraz ciąg operacji mnożenia tego wektora przez macierz A [ (^) n × n ]:
{ } { } { } { } { }
{ 1 } { } { } 0 1 { } 0
2 1 2 0
1 0
X A X A A ... A X A X
X A X A X
X A X
k (^) = ⋅ k = ⋅ ⋅ ⋅ ⋅ = k ⋅
M
Ponieważ każdy wektor w przestrzeni n wymiarowej można przedstawić jako kombinację liniową wektorów własnych naszej macierzy, to w szczególności:
=
n j j^ j
X α V 1
gdzie V j oznacza j −ty wektor własny a α j mnożnik. Przyjmijmy ponadto, dla prostoty zapisu, że
wartości własne są uporządkowane w kolejności od największej do najmniejszej co do modułu (wskaźnik j ). Podstawiając (58) do (57) 3 otrzymujemy:
{ } (^) ∑ ∑
=
n j k j j j
X k^ Ak V A V 1
1 1
Ponieważ, na mocy definicji problemu własnego, zachodzi związek:
Ak +^1 ⋅ V j = λ k j +^1 ⋅ V j , (60)
to ostatni człon wyrażenia (59) jest równy:
∑= ⋅^ +^ ⋅ = ∑= ⋅ +⋅
n j j
n j kj j j j Ak^ V V 1
1 1
co, po wyciągnięciu przed znak sumy dominującej wartości własnej ostatecznie prowadzi do nastę- pującej zależności, wiążącej wektor X { k^ +^1 }z wartościami i wektorami własnymi macierzy A :
{ } ⎥
∑=
k X k^ k j j V 1
1
1
1 11
WydziaPolitechnika Krakowskał Inżynierii Lądowej
Zgodnie z wcześniejszym założeniem λ 1 jest dominującą wartością własną, czyli zachodzi zależ-
ność:
1 dla 1 1
λ λ j < j ≠ , (63)
a więc gdy k →∞to 0
1 1
⎟⎟^ →
⎛ (^) j k +
λ . Jak z tego widać
{ } {^ { } k } skł
k sk^ kł k (^) X X V λ X^1 lim 1 , 1
→∞ =^ =.^ (64)
We wzorze (64) symbol skł oznacza dowolnie wybraną składową wektora. Ze względu na stabilność numeryczną procedury, po każdym kroku potęgowym (57) nowo uzyska- ny wektor X {^ k +^1 } normalizujemy, czyli zastępujemy wektorem W {^ i +^1 } o tym samym kierunku i zwrocie, ale długości równej 1, czyli W {^ i +^^1 }^ T^ ⋅ W { i^ +^1 }^ = 1 (gdyby tej operacji nie wykonać to długość wektora po kolejnych krokach albo rosła by do nieskończoności, jeżeli długość wektora startowego była większa od jedności, albo malała do zera w przeciwnym przypadku) a wartości własnej nie ob- liczamy wprost ze wzoru (64) ale z ilorazu Rayleigha:
{ } 1 { }{ }^ { }{ }= { } (^) ⋅ { + 1 } ⋅ = ⋅ ⋅ kT k kT k k kT k λ XX AXX W X. (65)
Oczywiście po każdej iteracji należy sprawdzić warunki zakończenia obliczeń w postaci błędu zbieżności wartości własnej i wektora własnego: { } { } { } v { k^ }^ { } k^ v
k^ λ
k k λ w W W ε
ε λ
w λ λ
= − <
1
1
1
. (66)
Ponieważ tempo zbieżności wartości własnej jest znacznie wyższe niż tempo zbieżności wektora
własnego, dobór wartościε λ i ε v powinien być wykonywany rozważnie, zwłaszcza dla dużych
macierzy. W szczególności należy się zastanowić, czy w dalszej analizie są nam potrzebne zarówno wartość jak i odpowiadający jej wektor własny, a jeżeli tak, to czy muszą być wyznaczone tak samo precyzyjnie. Przy sprawdzaniu warunku (66) 2 należy mieć na uwadze fakt, że jeżeli dominująca wartość własna jest ujemna, to wektory będące kolejnymi przybliżeniami odpowiadającego jej wektora własnego co iteracja będą zmieniały zwrot na przeciwny. Ostateczny algorytm metody wygląda następująco:
Wyznaczyć wektor startowy X { }^0 mając na względzie fakt, że jego dobór może mieć znaczący wpływ na liczbę wykonanych iteracji. Oczywiście, jak łatwo się domyślić, im wektor startowy jest bliższy poszukiwanemu wektorowi własnemu, tym mniej iteracji trzeba będzie wykonać;
1. Znormalizować wyznaczony wektor { }^
{ } { } kT { } k
k k X X
W X
2. Wykonać krok potęgowy X { k^^ +^1 }^ = A ⋅ W { } k ; 3. Obliczyć iloraz Rayleigha λ { k^^ +^1 }^ = W { } kT^ ⋅ X {^ k +^1 };
WydziaPolitechnika Krakowskał Inżynierii Lądowej
- Warunki zakończenia obliczeń: { } { } { }
v { }^ { }^ (^ )^ v
λ λ
w W W ε
ε λ
λ λ w
−
−
− − 1 1 ,^141277
0 , 500
0 , 500
0 , 500
0 , 500
0 , 814
0 , 465
0 , 349
0 , 000 1 0
0
1 0
czyli konieczne jest wykonanie następnej iteracji. Iteracja trzecia:
- Normalizacja:
{ } { } { }
{ } ⎥
−
− −
− 0 , 759
0 , 517
0 , 348
0 , 190
8 , 370
5 , 696
3 , 836
2 , 092 2 2 2
2 11 , 026944
1 X 1
X X
W T. (75)
- Krok potęgowy:
{ } { } ⎥⎥
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢⎢
⎢
⎣
⎡ −
−
−
− − −
−
− − −
−
⎥
9 , 003
5 , 334
3 , 500
1 , 255
0 , 759
0 , 517
0 , 348
0 , 190
4 1 5 7
1 2 2 5
2 6 2 1
3 2 1 4 X^3 AW^2. (76)
- Iloraz Rayleigha:
{ } { } { } [ ] 11 , 044677
9 , 003
5 , 334
3 , 500
1 , 255 (^2 230) , 190 0 , 3480 , 517 0 , 759 =− ⎥
− λ W T^ X − −. (77)
- Warunki zakończenia obliczeń: { } { } { }
v { }^ { }^ (^ )^ v
λ λ
w W W ε
ε λ
w λ λ
− −
0 , 814
0 , 465
0 , 349
0 , 000
0 , 759
0 , 517
0 , 348
0 , 190 2 1
2
2 1
czyli dalej konieczne jest wykonanie następnej iteracji. W toku dalszych obliczeń można się przekonać, że po ośmiu iteracjach: { } { } { }
v { }^ { }^ (^ )^ v
λ λ
w W W ε
ε λ
λ λ w
−
−
−
0 , 808
0 , 489
0 , 292
0 , 150
0 , 808
0 , 490
0 , 291
0 , 154 8 7
8
8 7
czyli został spełniony warunek (66) 1 (błąd zbieżności ze względu na wartość własną), a po 15 itera- cjach:
WydziaPolitechnika Krakowskał Inżynierii Lądowej
{ } { } { }
v {^ }^ {^ }^ (^ )^ v
λ λ
w W W ε
ε λ
λ λ w
−
− −
− 1 0 , 000057
0 , 809
0 , 489
0 , 288
0 , 154
0 , 809
0 , 489
0 , 288
0 , 154 15 14
15
15 14
. (80)
czyli został spełniony warunek (66)^2 (błąd zbieżności ze względu na wektor własny). Jak widać, uzyskanie takiej samej dokładności w przypadku wektora własnego, jak w przypadku wartości wła- snej wymaga wykonania około dwukrotnie większej liczby iteracji.
Efektywność (tempo zbieżności) metody potęgowej zależy od wartości ilorazu λ 2 (^) λ 1 , gdyż
{ } { }
1
2 (^1112)
(^11112 11112 )
1
⎥ =^ ≤
≤ ⋅ ⋅ + ⋅⎛^ ⋅
=
= =
n j j j
k k
n j j j
k n k j j j j j
k k j
k k
c V c V
V V c V c V
X X
gdzie c j = ( λj − 1 ) ⋅ αj , j = 1 , 2 ,..., n są stałymi współczynnikami rozkładu. Ze wzoru (81) wynika, że
tempo, z jakim wektor X zmierza do wektora własnego V 1 jest szacowane od góry przez wartość tego ilorazu. Najgorsza sytuacja wystąpi wtedy, gdy iloraz jest równy jeden, czyli wartości te różnią się jedynie znakiem. W kolejnych iteracjach wystąpią oscylacje pomiędzy dwoma przeciwnymi wartościami λ bez szans na spełnienie warunku zakończenia obliczeń. Najlepsza sytuacja nato- miast (a zatem również największe przyspieszenie obliczeń) nastąpi wtedy, gdy iloraz ten będzie miał największą możliwą dla danej macierzy wartość. Sytuacja taka pojawi się, gdy λ 2 (^) = λn. Moż- na do niej doprowadzić przez odpowiednie przesunięcie widma wartości własnych wyjściowej ma- cierzy o wartość
λopt = ( λ 2 + λn ) 2. (82)
Należy w tym celu najpierw możliwie precyzyjnie oszacować wartości własne λ 2 i λn (na przykład korzystając z ciągów Sturma) a następnie przekształcić macierz A [ (^) n × n ] odejmując od wszystkich elementów jej przekątnej głównej składnik λopt
A '^ = A − λopt ⋅ I , (83)
gdzie I (^) [ n × n ]oznacza macierz jednostkową, i dalsze obliczenia prowadzić na tak uzyskanej macierzy
A [^ ' n × n ] , dla której będzie zachodził warunek λ' 2 = λ' n. Po zakończeniu iteracji należy oczywiście
powrócić do wartości własnej macierzy wyjściowej, dodając do wyniku obliczeń metodą potęgową λ opt. Wyznaczony w trakcie tych obliczeń wektor własny nie wymaga żadnych modyfikacji, gdyż przesunięcie widma nie ma wpływu na jego składowe.
WydziaPolitechnika Krakowskał Inżynierii Lądowej
1. Znormalizować wyznaczony wektor { }^
{ } { } kT { } k
k k X X
W X
2. Wykonać krok potęgowy A ⋅ X { k^^ +^1 }^ = W { } k przez rozwiązanie układu liniowych równań algebra- icznych; 3. Obliczyć iloraz Rayleigha λ { k^^ +^1 }^ = W { } kT^ ⋅ X {^ k +^1 }; 4. Sprawdzić warunek (66). W zależności od wyniku sprawdzenia albo należy wrócić do 1., albo zakończyć pracę. Jeżeli metodę iteracji odwrotnej zastosujemy do macierzy po przesunięciu widma wartości wła- snych o liczbę d (83), uzyskamy zbieżność rozwiązania nie do najbliższej zeru wartości własnej macierzy wyjściowej a do wartości własnej najbliższej tej zadanej liczbie. W ten sposób można, stosując wariant metody potęgowej, wyznaczyć dowolną wartość własną macierzy. Uogólniony problem własny Uogólnionym problemem własnym nazywamy zadanie w postaci: A ⋅ V = λ ⋅ B ⋅ V. (86) Jeżeli wyznacznik macierzy B jest różny od zera zadanie takie można formalnie sprowadzić do standardowego problemu własnego przez obustronne pomnożenie przez macierz B −^1 , co prowadzi do:
B −^1 ⋅ A ⋅ V = λ ⋅ V. (87) Takie postępowanie ma jednak zasadniczą wadę spowodowaną numeryczną złożonością procedury (odwracanie macierzy), oraz możliwością utraty stabilności w trakcie obliczeń. Ponadto, jeżeli ma- cierz B jest pasmowa (częsty przypadek) to macierz odwrotna do niej najczęściej taka nie jest, co w przypadku dużych zadań może być kłopotliwe. Jeżeli macierz B jest symetryczna i dodatnio określona możliwy jest następujący sposób postępo- wania, polegający na zastosowaniu rozkładu L ⋅ LT :
B = L ⋅ L^ T , (88) co po podstawieniu do (86) daje:
A ⋅ X =λ ⋅ L ⋅ LT^ ⋅ X , (89)
a po lewostronnym pomnożeniu przez L −^1 :
L −^1 ⋅ A ⋅ X =λ ⋅ L −^1 ⋅ L ⋅ LT^ ⋅ X =λ⋅ LT ⋅ X. (90)
Jeżeli dokonamy teraz podstawienia:
LT ⋅ X = Y , czyli X = L − T ⋅ Y , (91)
to problem (86) sprowadzimy do standardowego problemu własnego (1), z wartościami własnymi zachowanymi z zadania (86), ale zmienionymi wektorami własnymi:
Y Y C
L 1 − 41 ⋅ 2 A ⋅ 4 L 3 − T^ ⋅ =λ⋅ , (92)
który możemy rozwiązać na przykład jedną z przedstawionych powyżej metod. Znając wektory własne problemu (92) możemy w każdej chwili wrócić do wektorów własnych zadania wyjściowe- go (86), korzystając z zależności (91). Jeżeli natomiast w trakcie rozwiązywania zadania (86) meto-
WydziaPolitechnika Krakowskał Inżynierii Lądowej
dą potęgową przy uwzględnieniu zależności (88) chcemy uniknąć konieczności odwracania macie- rzy L i LT możemy zastosować procedurę postępowania przedstawioną poniżej:
Dokonać rozkładu L ⋅ LT macierzy B i przyjąć wektor startowy Y {^0 };
1. Znormalizować wyznaczony wektor { }^
{ } { } kT { } k
k k Y Y
W Y
2. Obliczyć wektor X { } k rozwiązując układ liniowych równań algebraicznych LT ⋅ X { } k^ = W { } k ; 3. Obliczyć wektor Y {^ k +^1 } rozwiązując układ liniowych równań algebraicznych L ⋅ Y { k^^ +^1 }^ = A ⋅ X { } k , co jest równoznaczne wykonaniu kroku potęgowego;
4. Obliczyć iloraz Rayleigha λ{ k +^^1 }^ = W { } kT^ ⋅ Y {^ k +^1 };
5. Sprawdzić warunek (66). Dokładnie tak samo jak we wcześniej rozważanych przykładach, jeżeli warunek (66) został speł- niony, kończymy obliczenia, w przeciwnym razie wracamy do 1. W tym miejscu warto zwrócić uwagę na fakt, że w trakcie obliczeń (punkt 2.) w każdej iteracji jest znajdowany wektor własny oryginalnego problemu własnego (86), tak więc nie należy go obliczać ponownie na zakończenie pracy. Przedstawiony powyżej sposób postępowania zastosujemy do wyznaczenia dominującej wartości własnej uogólnionego problemu własnego opisanego symetrycznymi macierzami A [ (^) 4 × 4 ] i B [ (^) 4 × 4 ] , z
dokładnością ελ = ε v = 0 , 0001 :
−
− −
− −
− −
− − −
−
8 10 2 46
2 5 9 2
2 10 5 10
4 2 2 8
2 2 1 3
3 6 4 1
1 3 6 2
2 1 3 2 A , B. (93)
Operacje wstępne – rozkład macierzy B i wektor startowy Y {^0 }:
−
−
− − − −
−
− −
− −
− 0 0 0 5
0 0 2 1
0 3 2 2
2 1 1 4
4 2 1 5
1 2 2 0
1 3 0 0
2 0 0 0
8 10 2 46
2 5 9 2
2 10 5 10
4 2 2 8 B L L^ T , { } ⎥
−
1
1
1
1 Y^0. (94)
Iteracja pierwsza:
- Normalizacja:
{ } { } { }
{ } ⎥
− −
−
0 , 500
0 , 500
0 , 500
0 , 500
1
1
1
1 0 0 0
0 4
1 Y 1
Y Y
W T. (95)
- Układ równań:
{ } { } { } { } ⎥
−
− −
−
−
−
− − 0 , 100
0 , 200
0 , 033
0 , 167 0 0 , 500
0 , 500
0 , 500
0 , 500 0 0 0 0 5
0 0 2 1
0 3 2 2
2 1 1 4 LT^ X^0 W^0 X X. (96)
WydziaPolitechnika Krakowskał Inżynierii Lądowej
Iteracja trzecia:
- Normalizacja:
{ } { } { }
{ } ⎥
−
−
−
−
−
−
0 , 543
0 , 629
0 , 302
0 , 467
1 , 233
1 , 428
0 , 686
1 , 060 2 2 2
2 5 , 153832
1 Y 1
Y Y
W T. (104)
- Układ równań:
{ } { } { } { } ⎥
−
−
−
−
−
−
− − 0 , 109
0 , 369
0 , 073
0 , 671 2 0 , 543
0 , 629
0 , 302
0 , 467 2 0 0 0 5
0 0 2 1
0 3 2 2
2 1 1 4 LT^ X^2 W^2 X X. (105)
- Układ równań:
{ } { } { } { } ⎥
−
−
− −
−
− −
− − −
−
− −
−
− 1 , 279
1 , 374
0 , 747
1 , 080 3 0 , 109
0 , 369
0 , 073
0 , 671
2 2 1 3
3 6 4 1
1 3 6 2
2 1 3 2 3 4 2 1 5
1 2 2 0
1 3 0 0
2 0 0 0 L Y^3 A X^2 Y Y. (106)
- Iloraz Rayleigha:
{ } { } { } [ ] 2 , 288566
1 , 279
1 , 374
0 , 747
1 , 080 (^2 230) , 467 0 , 302 0 , 629 0 , 543 = ⎥
−
−
−
λ W T Y − −. (107)
- Warunki zakończenia obliczeń: { } { } { }
v { }^ { }^ v
λ λ
w W W ε
ε λ
λ λ w
−
−
−
−
−
− 0 , 138843
0 , 559
0 , 551
0 , 246
0 , 567
0 , 543
0 , 629
0 , 302
0 , 467 2 1
1
2 1
czyli konieczne jest wykonanie następnej iteracji. W toku dalszych obliczeń można się przekonać, że po siedmiu iteracjach: { } { } { }
v { }^ { }^ v
λ λ
w W W ε
ε λ
λ λ w
−
−
−
−
−
− 0 , 011241
0 , 554
0 , 617
0 , 313
0 , 464
0 , 556
0 , 609
0 , 319
0 , 467 7 6
7
7 6
czyli został spełniony warunek (66) 1 (błąd zbieżności ze względu na wartość własną), a po 22 itera- cjach: { } { } { }
{ } { } (^) v ,
,
,
,
,
,
,
, v
λ λ
w W W ε
ε λ
λ λ w
−
−
−
−
− 0 , 000097
0555
0612
0317
0466
0555
0612
0317
0466 22 21
22
22 21
. (110)
WydziaPolitechnika Krakowskał Inżynierii Lądowej
czyli został spełniony warunek (66)^2 (błąd zbieżności ze względu na wektor własny). Jak widać, uzyskanie takiej samej dokładności w przypadku wektora własnego, jak w przypadku wartości wła- snej w tym przypadku wymaga wykonania około trzykrotnie większej liczby iteracji.
0 2 4 6 8 10 12 14 16 18 20 22
10 9
10 8
10 7
10 6
10 5
10 4
10 3
10
λn w (^) n
n n,
Rys. 3 Zbieżność metody potęgowej dla dominującej wartości własnej ( λ ) i odpowiadającego jej wektora własnego ( w ) w kolejnych iteracjach ( n ) uogólnionego problemu własnego, skala logarytmiczna.
