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:
V V
A⋅ =λ⋅ . (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 V1,V ...,2, V . Wektory te oczywiście nie n 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:
( ) ( )
0det A−λ⋅I = p λ = , (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:
( )
A =0p . (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 QT⋅Q=I to wartości własne macierzy QT⋅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:
max
0 λmin ≤λ
⋅
⋅
≤ ⋅
< 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:
• 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
[ ]n n
A × 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 .
WT
D W
A= ⋅ ⋅ , (5)
gdzie:
{ } ... Q{ } ... Q{ }2 Q{ }1
Q
W = m ⋅ ⋅ i ⋅ ⋅ ⋅ . (6)
jest iloczynem macierzy Q definiujących kolejne transformacje ortogonalne { }i 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:
1. 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:
{ } { }i
kl l
k n
l k i
pq a
a q p
≠
= = ,..., 2 , 1 ,max :
, . (7)
2. Wyznaczenie macierzy transformacji Q tak, aby elementy nowej macierzy: { }i
{ }i Q{ }iT A{ }i Q{ }i
A +1 = ⋅ ⋅ (8)
spełniały warunek:
{ }pqi+1 =aqp{ }i+1 =0
a , (9)
i wyznaczenie elementów nowej macierzy A{ }i+1 .Jak łatwo sprawdzić, jeżeli macierz Q { }i
s q
s q
c q q
qp pq
qq pp
−
=
=
=
=
, (10)
czyli ostateczna macierz będzie miała postać:
{ }
q p
q p i
c s
s Q c
⎥⎥
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎢⎢
⎣
⎡
= −
0 0 0 0
1 0 0
1 0
0
1 0 0 0
1
L L
L M M M
M L M L M L M M L M L M L M
M L L L
, (11)
to spełnienie warunku ortogonalności QT⋅Q=I jest równoznaczne z zażądaniem aby wielkości oznaczone jako c i s spełniały zależność:
2 1
2+ s =
c . (12)
Warunek ten jest jednocześnie warunkiem wystarczającym ortogonalności macierzy Q .
Jak łatwo się przekonać, w trakcie transformacji (8) zmieniają się jedynie elementy macierzy A należące do wierszy i kolumn opatrzonych wskaźnikami p i q . Odpowiednie wzory transfor- macyjne mają postać:
{ } { } { } { }
{ } { } { } { }
{ }pqi
( )
{ }pqi(
{ }ppi { }qqi)
i pq i
qq i
pp i
i pq i
qq i
pp i
pp
a a s c a s c a
a c s a
c a s a
a c s a
s a c a
−
⋅
⋅ +
⋅
−
=
⋅
⋅
⋅ +
⋅ +
⋅
=
⋅
⋅
⋅
−
⋅ +
⋅
=
+ + +
2 2 1
2 2
1
2 2
1
2 2
, (13)
{ } { } { }
{ } { } { }
{ } { } { }
{ } { } { }
q r p r a c a s a
a s a c a
a c a s a
a s a c a
i qr i
pr i
rq
i qr i
pr i
rq
i rq i
rp i
rq
i rq i
rp i
rp
≠
≠
⋅ +
⋅
=
⋅
−
⋅
=
⋅ +
⋅
=
⋅
−
⋅
=
+ + + +
,
1 1 1 1
, (14)
i można je wyprowadzić na przykład przez bezpośrednie wykonanie operacji macierzowych wy- stępujących w (8).
Warunek (9) w połączeniu z równaniem (13)3 służy do wyznaczenia c i s dla bieżącej iteracji.
A mianowicie, ponieważ:
(
2− 2)
⋅ { }+ ⋅ ⋅(
{ }− qq{ }i)
=0i pp i
pq c s a a
a s
c (15)
jest równoważne:
{ } { } { }i pq i pp i qq
a a a s c
s c
⋅
= −
⋅
⋅
= −
2 2
2
η 2 , (16)
oznaczonemu η , przy podstawieniu
c
t = s (17)
równanie (16) przechodzi w:
0 1
2+2⋅ ⋅t− =
t η , (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łę:
( )
1 sgn
2+
= +
η η
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łę:
η
= ⋅ 2
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):
c t s t c
⋅
= +
= 1
1
2 . (21)
3. Obliczenie macierzy wektorów własnych na podstawie zależności:
{ }i W{ }i Q{ }i
W +1 = ⋅ , (22)
przy czym początkową macierzą W jest macierz jednostkowa. Stosowne wzory transforma-{ }0 cyjne przyjmą postać:
{ } { } { }
{ } { } { } r n
w c w s w
w s w c w
i rq i
rp i
rq
i rq i
rp i
rp1 1,2,...,
1
⋅ = +
⋅
=
⋅
−
⋅
=
+ +
, (23)
która wynika z faktu, że w kolejnych iteracjach zmieniają się wyłącznie kolumny p i q macie- rzy W .
4. Sprawdzenie warunku zakończenia obliczeń:
{ } { }ii++11 <ε sn
sp , (24)
w którym licznik i mianownik ułamka są odpowiednio równe:
{ } { }
{ } { }1
,..., 2 , 1 1
1 ,..., 2 , 1 , 1
max max
+
= +
+
≠
= +
=
=
i n ii i i
i ij j
i n
j i i
a sn
a sp
(25)
dla normy maksimum, lub:
{ }
( )
{ }{ }
∑ ( ){ }
∑∑
+ +
=≠ = +
− +
=
=
n i
i
j i n i
n j
i n ij
n i
a sn
a sp
1 2 1 1
1 1
1 2 1 1
2
(26)
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 leżących poza przekątną główną i oznaczymy { }i ją przez S : { }i
{ }
( )
{ }j i n i
n j
i ij
i a
S
=≠ =
∑∑
=
1 1
2 , (27)
to wprost z (14) wynika, że:
{ }i 1 S{ }i 2
( )
a{ }pqi 2S + = − ⋅ , (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⋅ do wyznaczenia wartości i wektorów własnych tak zredukowanej R 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ą
0001 ,
=0
ε :
{ }
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
=
−
−
−
−
−
−
−
7 5 1 4
5 2 2 1
1 2 6 2
4 1 2 3
A0
A . (29)
Iteracja pierwsza:
1. Element wiodący (największy co do modułu element macierzy A spoza przekątnej głównej) { }0
=3
p , 4q= .
2. Współczynniki η i t :
{ } { } { }
{ } { } { }
(
0,500000)
1 0,618034 500000, 0
) 500000 ,
0 sgn(
1 ) sgn(
500000 ,
5 0 2
2 7 2
7 ,
5 , 2
2 2
0 34 0 33 0 44
0 44 0
34 0
33
− + =
− +
−
= − +
= +
−
⋅ = +
= −
⋅
= −
−
=
=
−
=
η η t η
a a η a
a a
a
. (30)
3. Współczynniki c i s :
( )
525731 ,
0 850651 ,
0 618034 ,
0
850651 ,
1 0 618034 ,
0 1 1
1
2 2
−
=
⋅
−
=
⋅
=
+ =
= −
= + c t s c t
. (31)
4. Macierz A : { }1
{ }
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
−
= −
851 , 0 526 , 0 000 , 0 000 , 0
526 , 0 851 , 0 000 , 0 000 , 0
000 , 0 000 , 0 000 , 1 000 , 0
000 , 0 000 , 0 000 , 0 000 , 1 0
0 0
0 0
0 0 1 0
0 0 0 1
c s
s
Q c . (32)
{ }1 =Q{ }0 ⋅A{ }0 ⋅Q{ }0 =
A T
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
⋅⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
⋅⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
090 , 10 000 , 0 902 , 1 877 , 2
000 , 0 090 , 1 176 , 1 954 , 2
902 , 1 176 , 1 000 , 6 000 , 2
877 , 2 954 , 2 000 , 2 000 , 3
526 , 0 526 , 0 000 , 0 000 , 0
851 , 0 851 , 0 000 , 0 000 , 0
000 , 0 000 , 0 000 , 1 000 , 0
000 , 0 000 , 0 000 , 0 000 , 1
2 5 1 4
5 7 2 1
1 2 6 2
4 1 2 3
851 , 0 526 , 0 000 , 0 000 , 0
526 , 0 851 , 0 000 , 0 000 , 0
000 , 0 000 , 0 000 , 1 000 , 0
000 , 0 000 , 0 000 , 0 000 , 1
. (33)
5. Macierz W : { }1
{ }1 =W{ }0 ⋅Q{ }0 = W
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
⋅⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡ − −
851 , 0 526 , 0 000 , 0 000 , 0
526 , 0 851 , 0 000 , 0 000 , 0
000 , 0 000 , 0 000 , 1 000 , 0
000 , 0 000 , 0 000 , 0 000 , 1
851 , 0 526 , 0 000 , 0 000 , 0
526 , 0 851 , 0 000 , 0 000 , 0
000 , 0 000 , 0 000 , 1 000 , 0
000 , 0 000 , 0 000 , 0 000 , 1
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
. (34)
6. Warunek zakończenia obliczeń (w normie maksimum):
{ } { }
{ } max { } 10,090170 953575 ,
2 max
1 4 , 3 , 2 , 1 1
1 4 , 3 , 2 , 1 , 1
=
=
=
=
=
= ≠
i ii
j i j ij
i
a sn
a sp
, (35)
{ }
{ } ε
sn
sp = =0,292718>
090170 ,
10
953575 ,
2
1 1
, (36)
czyli konieczne jest wykonanie następnej iteracji.
Iteracja druga:
1. Element wiodący (największy co do modułu element macierzy A spoza przekątnej głównej) { }1
=1
p , 3q= .
2. Współczynniki η i t :
{ } { } { }
{ } { } { }
(
0,323308)
1 0,727657 323308, 0
) 323308 ,
0 sgn(
1 ) sgn(
323308 ,
953575 0 ,
2 2
3 090170 ,
1 2
090170 ,
1 , 953575 ,
2 , 3
2 2
1 13 1 11 1 33
1 33 1
13 1
11
− + =
− +
−
= − +
= +
−
⋅ =
= −
⋅
= −
=
=
=
η η t η
a a η a
a a
a
. (37)
3. Współczynniki c i s :
( )
588375 ,
0 808588 ,
0 727657 ,
0
808588 ,
1 0 727657 ,
0 1 1
1
2 2
−
=
⋅
−
=
⋅
=
+ =
= −
= + c t s c t
. (38)
4. Macierz A{ }2 :
{ }
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
= −
−
000 , 1 000 , 0 000 , 0 0
000 , 0 809 , 0 000 , 0 588 , 0
000 , 0 000 , 0 000 , 1 000 , 0
000 , 0 588 , 0 000 , 0 809 , 0 1
1 0 0 0
0 0
0 0 1 0
0 0
c s
s c
Q . (39)
{ }2 =Q{ }1 ⋅A{ }1 ⋅Q{ }1 =
A T
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
⋅⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
⋅⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
090 , 10 693 , 1 902 , 1 326 , 2
693 , 1 059 , 1 127 , 2 000 , 0
902 , 1 127 , 2 000 , 6 926 , 0
326 , 2 000 , 0 926 , 0 149 , 5
000 , 1 000 , 0 000 , 0 000 , 0
000 , 0 809 , 0 000 , 0 588 , 0
000 , 0 000 , 0 000 , 1 000 , 0
000 , 0 588 , 0 000 , 0 809 , 0
090 , 10 000 , 0 902 , 1 877 , 2
000 , 0 090 , 1 176 , 1 954 , 2
902 , 1 176 , 1 000 , 6 000 , 2
877 , 2 954 , 2 000 , 2 000 , 3
000 , 1 000 , 0 000 , 0 000 , 0
000 , 0 809 , 0 000 , 0 588 , 0
000 , 0 000 , 0 000 , 1 000 , 0
000 , 0 588 , 0 000 , 0 809 , 0
. (40)
5. Macierz W{ }2 :
{ }2 =W{ }1 ⋅Q{ }1 = W
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
⋅⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡ −
−
−
−
851 , 0 425 , 0 000 , 0 309 , 0
526 , 0 688 , 0 000 , 0 501 , 0
000 , 0 000 , 0 000 , 1 000 , 0
000 , 0 588 , 0 000 , 0 809 , 0
000 , 1 000 , 0 000 , 0 000 , 0
000 , 0 809 , 0 000 , 0 588 , 0
000 , 0 000 , 0 000 , 1 000 , 0
000 , 0 588 , 0 000 , 0 809 , 0
851 , 0 526 , 0 000 , 0 000 , 0
526 , 0 851 , 0 000 , 0 000 , 0
000 , 0 000 , 0 000 , 1 000 , 0
000 , 0 000 , 0 000 , 0 000 , 1
. (41)
6. Warunek zakończenia obliczeń (w normie maksimum):
{ } { }
{ } max { } 10,090170 326205 ,
2 max
2 4 , 3 , 2 , 1 2
2 4 , 3 , 2 , 1 , 2
=
=
=
=
=
=≠
i ii j i j ij i
a sn
a sp
, (42)
{ }
{ } ε
sn
sp = =0,230542>
090170 ,
10
326205 ,
2
2 2
, (43)
czyli dalej konieczne jest wykonanie następnej iteracji.
Iteracja trzecia:
1. Element wiodący (największy co do modułu element macierzy A spoza przekątnej głównej) { }2
=1
p , 4q= .
2. Współczynniki η i t :
{ } { } { }
{ } { } { }
(
3,275584)
1 0,149245 275584, 3
) 275584 ,
3 sgn(
1 ) sgn(
275584 ,
326205 3 ,
2 2
149190 ,
5 090170 ,
10 2
090170 ,
10 ,
326205 ,
2 ,
149190 ,
5
2 2
2 14 2 11 2 44
2 44 2
14 2
11
− + =
− +
−
= − +
= +
−
⋅ =
−
= −
⋅
= −
−
=
=
=
η η t η
a a η a
a a
a
. (44)
3. Współczynniki c i s :
( )
147610 ,
0 989046 ,
0 149245 ,
0
989046 ,
1 0 149245 ,
0 1 1
1
2 2
−
=
⋅
−
=
⋅
=
+ =
= −
= + c t s c t
. (45)
4. Macierz A{ }3 :
{ }
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
−
=
−
989 , 0 000 , 0 000 , 0 148 , 0
000 , 0 000 , 0 000 , 0 000 , 0
000 , 0 000 , 0 000 , 1 000 , 0
148 , 0 000 , 0 000 , 0 989 , 0 2
0 0
0 1 0 0
0 0 1 0
0 0
c s
s c
Q . (46)
{ }3 =Q{ }2 ⋅A{ }2 ⋅Q{ }2 =
A T
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
⋅⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
⋅⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
437 , 10 674 , 1 745 , 1 000 , 0
674 , 1 059 , 1 127 , 2 250 , 0
745 , 1 127 , 2 000 , 6 196 , 1
000 , 0 250 , 0 196 , 1 496 , 5
989 , 0 000 , 0 000 , 0 148 , 0
000 , 0 000 , 1 000 , 0 000 , 0
000 , 0 000 , 0 000 , 1 000 , 0
148 , 0 000 , 0 000 , 0 989 , 0
090 , 10 693 , 1 902 , 1 326 , 2
693 , 1 059 , 1 127 , 2 000 , 0
902 , 1 127 , 2 000 , 6 926 , 0
326 , 2 000 , 0 926 , 0 149 , 5
989 , 0 000 , 0 000 , 0 148 , 0
000 , 0 000 , 1 000 , 0 000 , 0
000 , 0 000 , 0 000 , 1 000 , 0
148 , 0 000 , 0 000 , 0 989 , 0
. (47)
5. Macierz W : { }3
{ }3 =W{ }2 ⋅Q{ }2 = W
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
⋅⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡ −
−
−
−
−
−
796 , 0 425 , 0 000 , 0 432 , 0
594 , 0 688 , 0 000 , 0 417 , 0
000 , 0 000 , 0 000 , 1 000 , 0
119 , 0 588 , 0 000 , 0 800 , 0
989 , 0 000 , 0 000 , 0 148 , 0
000 , 0 000 , 1 000 , 0 000 , 0
000 , 0 000 , 0 000 , 1 000 , 0
148 , 0 000 , 0 000 , 0 989 , 0
851 , 0 425 , 0 000 , 0 309 , 0
526 , 0 688 , 0 000 , 0 501 , 0
000 , 0 000 , 0 000 , 1 000 , 0
000 , 0 588 , 0 000 , 0 809 , 0
. (48)
6. Warunek zakończenia obliczeń (w normie maksimum):
{ } { }
{ } max { } 10,437343 127302 ,
2 max
3 4 , 3 , 2 , 1 3
3 4 , 3 , 2 , 1 , 3
=
=
=
=
=
=≠
i ii j i j ij i
a sn
a sp
, (49)
{ }
{ } ε
sn
sp = =0,203816>
437343 ,
10
127302 ,
2
3 3
, (50)
czyli dalej konieczne jest wykonanie następnej iteracji.
Ostatecznie po wykonaniu 11 iteracji otrzymujemy:
{ }
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
≅
−
−
−
−
137 , 11 000 , 0 001 , 0 000 , 0
000 , 0 103 , 0 000 , 0 000 , 0
001 , 0 000 , 0 626 , 6 000 , 0
000 , 0 000 , 0 000 , 0 661 , 5
A11
D , (51)
{ }
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
≅
−
−
−
−
−
−
809 , 0 276 , 0 318 , 0 410 , 0
489 , 0 794 , 0 109 , 0 346 , 0
288 , 0 350 , 0 883 , 0 119 , 0
154 , 0 414 , 0 327 , 0 836 , 0
W 11
W , (52)
i warunek zakończenia obliczeń:
{ } { }
{ } max { } 11,137200 000265 ,
0 max
11 4 , 3 , 2 , 1 11
11 4 , 3 , 2 , 1 , 11
=
=
=
=
=
=≠
i ii j i j ij i
a sn
a sp
, (53)
{ }
{ } ε
sn
sp = =0,000067<
137200 ,
11
000265 ,
0
11 11
, (54)
podczas gdy wartości ścisłe macierzy D i W są równe odpowiednio:
{ }
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
=
−
−
137 , 11 000 , 0 000 , 0 000 , 0
000 , 0 103 , 0 000 , 0 000 , 0
000 , 0 000 , 0 626 , 6 000 , 0
000 , 0 000 , 0 000 , 0 661 , 5
Aść
D , (55)
{ }
⎥⎥
⎦
⎤
⎢⎢
⎣
=⎡
=
−
−
−
−
−
−
−
−
−
−
809 , 0 276 , 0 318 , 0 410 , 0
489 , 0 794 , 0 109 , 0 346 , 0
288 , 0 350 , 0 883 , 0 119 , 0
154 , 0 414 , 0 327 , 0 836 , 0
W ść
W . (56)
Jak widać przyjęcie warunku zakończenia obliczeń w normie maksimum z dopuszczalnym pozio- mem błędu 0001ε =0, pozwoliło na wyznaczenie wszystkich wartości i wektorów własnych ma- cierzy A z dokładnością co najmniej trzech cyfr znaczących.
0 2 4 6 8 10 12
0.3 0.2 0.1 0.1 0.2 0.3
zk
k
Rys. 1. Błąd metody Jacobiego w normie maksimum (k – liczba iteracji).
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 oraz ciąg operacji mnożenia tego wektora przez macierz { }0 A[ ]n×n :
{ } { }
{ } { } { }
{ }1 { } { }0 1 { }0
0 2 1 2
0 1
... A X A X
A A X A X
X A X A X
X A X
k k
k = ⋅ = ⋅ ⋅ ⋅ ⋅ = ⋅
⋅
=
⋅
=
⋅
=
+
+ M . (57)
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 Vj
X
1
0 , (58)
gdzie Vj 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 j
k j n
j j j
k
k A V A V
X
1
1 1
1
1 α α . (59)
Ponieważ, na mocy definicji problemu własnego, zachodzi związek:
j k j j
k V V
A +1⋅ =λ+1⋅ , (60)
to ostatni człon wyrażenia (59) jest równy:
∑
∑
=+
=
+ ⋅ = ⋅ ⋅
⋅ n
j j
k j j n
j j
k
j A V V
1
1 1
1 α λ
α , (61)
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 :
{ }
⎥⎥
⎦
⎤
⎢⎢
⎣
⎡ ⎟⎟⎠ ⋅
⎜⎜ ⎞
⎝
⋅⎛
⋅
=
∑
=
+ +
+ n
j j
k j j k
k V
X
1
1
1 1
1 1
λ α λ
λ . (62)
