Numeryczne zadanie własne polega na znalezieniu niektórych lub wszystkich wartości oraz wektorów własnych macierzy wymiaru
rzeczywistej lub zespolonej. Na wykładzie z Matematyki Obliczeniowej omawia się kilka najprostszych metod dla symetrycznego zadania własnego, tzn. dla zadania własnego z macierzą rzeczywistą symetryczną.
Funkcją octave'a służącą do znajdowania wartości i wektorów własnych dowolnej macierzy jest
eig()
.
Najprostsze jej wywołanie to:
po którym funkcja powinna zwrócić wszystkie wartości własne w wektorze .
Jeśli potrzebujemy poznać wektory własne, to musimy wywołać tę funkcję z dwoma zwracanymi wartościami:
W tym przypadku będzie macierzą diagonalną, na której diagonali będą wartości własne macierzy
, natomiast
będzie macierzą ortogonalną, w której kolumnach będą wektory własne dla odpowiednich wartości na diagonali
, czyli
Sprawdźmy na losowej macierzy diagonalnej jak działa funkcja eig()
,
najpierw przy prostszym wywołaniu:
Wyznacznik jest na poziomie
, czyli de facto jest równy zero.
Teraz sprawdźmy drugą formę wywołania funkcji eig()
:
Sprawdźmy, jak działa funkcja dla macierzy większego wymiaru:
Wniosek: funkcja działa dobrze nawet dla większych .
Przetestujmy ją dla macierzy ogromnego wymiaru o dużej ilości zer, np. macierzy trójdiagonalnej. Wykorzystamy też funkcje sparse()
służące utworzeniu macierzy w formacie rzadkim:
Sprawdźmy, jak działa funkcja dla macierzy niesymetrycznej:
Istnieje też funkcja eigs()
, która oblicza tylko kilka wartości własnych macierzy, np. największych co do modułu wartości.
Przetestujmy funkcję dla losowej macierzy wymiaru :
Funkcja domyślnie znajduje sześć największych co do modułu wartości własnych, ale może znaleźć mniej lub więcej:
lub bazując na innym kryterium, np. najmniejsze wartości własne co do modułu:
Możliwe opcje to:
Spróbujmy policzyć dwie największe i dwie najmniejsze co do wartości własne macierzy trójdiagonalnej przy pomocy tej funkcji:
Niestety metoda zawiodła, a funkcja eig()
- choć wolna, zwróciła wszystkie wartości i wektory tej macierzy.
Zaimplementujmy najprostszą wersję metody potęgowej:
Niech będzie wektorem o normie drugiej równej jeden. Wtedy
ciąg iteracji metody potęgowej jest zdefiniowany następująco:
\begin{equation}(#)
\begin{array}{lclr}
\hat{x}_k&=&Ax_{k-1} &\qquad \mathrm{iteracja}\\
x_k&=&\frac{ \hat{x}_k}{ \|\hat{x}_k\|_2} & \qquad \mathrm{normowanie} \\
r_k&=& x_k^TAx_k \qquad &\mathrm{iloraz Rayleigha}
\end{array} \qquad k>0
\end{equation}
W octavie kod może wyglądać następująco:
Za warunek stopu przyjęliśmy zadanej tolerancji.
Przetestujmy ten algorytm dla macierzy symetrycznej A=[4,1;1,4]
, która ma wartości własne .
Dla TOL=1e-12
otrzymaliśmy, że po iteracjach iloraz Rayleigha wynosi
, czyli błąd
.
Z kolei przybliżenie wektora własnego to
0.707106628756140
0.707106933616922
Możemy policzyć w octave'ie, czyli
norm(A*x-5*x,2)
, która wynosi
.
Po prostych rachunkach otrzymujemy, że wektorem własnym dla wartości własnej jest wektor postaci
.
Sprawdźmy, dla naszego wektora:
Oczywiście ilość iteracji i błędy zależą też od wektora .
Powtórzmy obliczenia drukując na ekranie błąd i
dla
iteracji metody potęgowej i
ilorazu Rayleigha:
4 1
1 4
[it]| Rayleigh | |x-y| |
-----------------------------
[ 1] | 2.5008722 | 0.0241133 |
[ 2] | 4.9995814 | 0.0204587 |
[ 3] | 4.9998493 | 0.0122760 |
[ 4] | 4.9999457 | 0.0073658 |
[ 5] | 4.9999805 | 0.0044195 |
[ 6] | 4.9999930 | 0.0026517 |
[ 7] | 4.9999975 | 0.0015910 |
[ 8] | 4.9999991 | 0.0009546 |
[ 9] | 4.9999997 | 0.0005728 |
[10] | 4.9999999 | 0.0003437 |
[11] | 5.0000000 | 0.0002062 |
[12] | 5.0000000 | 0.0001237 |
[13] | 5.0000000 | 0.0000742 |
[14] | 5.0000000 | 0.0000445 |
[15] | 5.0000000 | 0.0000267 |
[16] | 5.0000000 | 0.0000160 |
[17] | 5.0000000 | 0.0000096 |
[18] | 5.0000000 | 0.0000058 |
[19] | 5.0000000 | 0.0000035 |
[20] | 5.0000000 | 0.0000021 |
[21] | 5.0000000 | 0.0000012 |
-----------------------------
x =
7.07106556762792e-01
7.07107005610232e-01
Wartosc wlasna: 4.999999999999440
Powtarzając eksperyment dla macierzy , czyli macierzy o wartościach własnych
i
, otrzymujemy po
iteracjach iloraz Rayleigha taki, że błąd
wynosi
przy tym samym warunku stopu:
2 1
1 2
[it]| Rayleigh+p | |x-y| |
-----------------------------
[ 1] | 3.5005623 | 0.0335338 |
[ 2] | 4.9997501 | 0.0158070 |
[ 3] | 4.9999722 | 0.0052693 |
[ 4] | 4.9999969 | 0.0017564 |
[ 5] | 4.9999997 | 0.0005855 |
[ 6] | 5.0000000 | 0.0001952 |
[ 7] | 5.0000000 | 0.0000651 |
[ 8] | 5.0000000 | 0.0000217 |
[ 9] | 5.0000000 | 0.0000072 |
[10] | 5.0000000 | 0.0000024 |
[11] | 5.0000000 | 0.0000008 |
-----------------------------
x =
7.07106736568327e-01
7.07106825804765e-01
Wartosc wlasna: 4.999999999999928
Widać już z tych dwóch eksperymentów, że ilość iteracji zależy od stosunku ( gdzie
to wartość własna o maksymalnym module,
to druga z kolei wartość własna ze względu na moduł) i że szybkość zbieżności ilorazu Rayleigha jest dużo większa niż
wektorów iteracji do odpowiedniego wektora własnego.
Potwierdza to oszacowania teoretyczne, w których można pokazać, że
Widać, że iloraz Rayleigha zbiega do wartości własnej dużo szybciej w tym przypadku.
Jeśli to para własna dla macierzy
, to
jest parą własną dla macierzy
. Ta obserwacja pozwala nam skonstruować tzw. odwrotną metodę potęgową, która jest zdefiniowana jako metoda potęgowa zastosowana do macierzy
dla
parametru będącego dobrym przybliżeniem
dowolnej wartości własnej
. Warunkiem zbieżności iteracji odwrotnej metody potęgowej do
wektora własnego dla
, przy założeniu, że wektor startowy nie jest ortogonalny do
jest to, aby
Szybkość zbieżności wektorów iteracji dla odwrotnej metody potęgowej zależy wtedy od
.
W octave'ie łatwo zaimplementować odwrotną metodę potęgową: wystarczy mnożenie przez macierz w metodzie potęgowej zastąpić przez mnożenie wektora przez macierz odwrotną, które jest równoważne rozwiązywaniu układu równań. W octave kod odwrotnej metody potęgowej może wyglądać następująco:
Zobaczmy jak szybko działa metoda dla naszej testowej macierzy i
:
-1.0100 1.0000
1.0000 -1.0100
[it]| 1/r+p | |x-y| |
-----------------------------
[ 1] | 4.9900001 | 0.0273613 |
[ 2] | 5.0000000 | 0.0001925 |
[ 3] | 5.0000000 | 0.0000010 |
[ 4] | 5.0000000 | 0.0000000 |
-----------------------------
x =
-7.07106781186489e-01
-7.07106781186606e-01
Wartosc wlasna: 5.000000000000000
A teraz dla - tu wektor własny ma postać
:
1.2000 1.0000
1.0000 1.2000
[it]| 1/r+p | |x+y| |
-----------------------------
[ 1] | 7.1576796 | 1.0000000 |
[ 2] | 3.9194920 | 1.3442338 |
[ 3] | 3.0139676 | 0.3789639 |
[ 4] | 3.0001162 | 0.0357476 |
[ 5] | 3.0000010 | 0.0032508 |
[ 6] | 3.0000000 | 0.0002955 |
[ 7] | 3.0000000 | 0.0000269 |
[ 8] | 3.0000000 | 0.0000024 |
[ 9] | 3.0000000 | 0.0000002 |
-----------------------------
x =
7.07106782104050e-01
-7.07106780269045e-01
Wartosc wlasna: 3.000000000000000
Przesuwanie widma macierzy poprzez dodanie do niej jest bardzo ważne w różnych uogólnieniach metody potęgowej takich jak np. metoda
z przesunięciami, por. np. [Kincaid:2006:AN], [Trefethen:1997:NLA], [Jankowska:1982:MN].