|
|
(Nie pokazano 19 wersji utworzonych przez 4 użytkowników) |
Linia 1: |
Linia 1: |
| ==Wyznaczanie wektorów i wartości własnych== | | == Flash == |
| | <flash>file=Test19.swf|width=750|height=530</flash> |
| | http://osilek.mimuw.edu.pl/images/2/29/Test19.swf |
|
| |
|
| {{przyklad|Moj przykład||Tekst}}
| | ===Metoda iteracji prostej Banacha=== |
| {{algorytm|Nie robiący nic||Leż}}
| |  |
| | Zupełnie inne, i jak się okaże --- przy odrobinie sprytu bardzo skuteczne --- |
| | podejście do wyznaczania miejsca zerowego jest oparte na metodzie Banacha. |
|
| |
|
| Niech będzie dana rzeczywista kwadratowa macierz <math>\displaystyle A</math> wymiaru <math>\displaystyle N</math>. Wektorem własnym <math>\displaystyle x\in C^N</math> oraz
| | Najpierw nasze równanie nieliniowe |
| odpowiadającą mu wartością własną <math>\displaystyle \lambda \in C</math> nazwiemy taką parę, dla której
| |
|
| |
|
| <center><math>\displaystyle A x = \lambda x, | | <center><math> |
| | f(x) = 0 |
| </math></center> | | </math></center> |
|
| |
|
| przy czym <math>\displaystyle x\neq 0</math>.
| | przekształcamy (dobierając odpowiednią funkcję <math>\phi</math>) do równania równoważnego |
| | (tzn. mającego te same rozwiązania) |
|
| |
|
| Zadanie wyznaczania wartości własnych i wektorów własnych macierzy ma bardzo
| | <center><math> |
| szerokie zastosowania w tak odległych do siebie dziedzinach jak np. analiza
| | x\,=\,\phi( x). |
| odporności konstrukcji mechanicznych (wieżowce, mosty, wagony kolejowe) na
| |
| wibracje, czy też rankingowanie stron internetowych w wyszukiwarce Google.
| |
| | |
| {{przyklad|Odporność budynku na trzęsienie ziemi||
| |
| | |
| Rozważmy prosty układ mechaniczny opisujący, naturalnie w pewnym
| |
| jedynie przybliżeniu, zachowanie się układu <math>\displaystyle N</math> ciężkich płyt połączonych ze
| |
| sobą relatywnie elatycznymi dźwigarami --- co może np. modelować konstrukcję
| |
| wieżowca.
| |
| | |
| {}{Model wieżowca poddanego drganiom poprzecznym}
| |
| | |
| Wiadomo, że jeśli częstotliwości drgań własnych tego wieżowca będą bliskie
| |
| częstotliwości siły wymuszającej (o niewielkiej amplitudzie), to konstrukcja
| |
| wpadnie w rezonans i w końcu rozpadnie się wskutek zbyt wielkich przemieszczeń.
| |
| Wychylenia naszych płyt z położenia równowagi są opisywane układem pewnych
| |
| równań różniczkowych.
| |
| Teoria matematyczna takich równań różniczkowych pokazuje, że częstotliwości
| |
| drgań własnych to nic innego jak ''wartości własne'' pewnej
| |
| {niesymetrycznej} macierzy
| |
| wymiaru <math>\displaystyle 2N</math>,
| |
| która powstaje ze współczynników równania różniczkowego opisującego dynamikę
| |
| tego układu.
| |
| }}
| |
| | |
| {{przyklad|Macierz Google'a||
| |
| | |
| Podstawowy algorytm rankingowania stron WWW w {http://www.wikipedia.org/pagerank}{wyszukiwarce Google}
| |
| sprowadza się do znalezienia rzeczywistego ''wektora własnego'' <math>\displaystyle \pi</math> pewnej silnie
| |
| rozrzedzonej macierzy <math>\displaystyle A</math> (gigantycznego rozmiaru, równego liczbie indeksowanych
| |
| stron, czyli w chwili pisania tego tekstu około {ilu} stron), odpowiadającego wartości własnej równej 1:
| |
| | |
| <center><math>\displaystyle A \pi = \pi.
| |
| </math></center> | | </math></center> |
|
| |
|
| Współrzędne wektora <math>\displaystyle \pi</math>
| | Następnie, startując z pewnego przybliżenia |
| interpretuje się jako wartość rankingową kolejnych stron WWW. Aby wszystko miało
| | początkowego <math>x_0</math>, konstruujemy ciąg kolejnych |
| sens, współrzędne wektora muszą być z przedziału [0,1]. Pewne
| | przybliżeń <math>x_k</math> według wzoru |
| twierdzenia matematyczne i subtelny dobór macierzy <math>\displaystyle A</math> gwarantują, że taki
| |
| wektor <math>\displaystyle \pi</math> zawsze istnieje i jest jedyny! Co więcej, wartość 1 jest
| |
| dominującą wartością własną <math>\displaystyle A</math>, a to z kolei ma ważne znaczenie dla tzw.
| |
| {sec:metoda-potegowa}{metody potęgowej} numerycznego wyznaczania takiego wektora.
| |
| }}
| |
| | |
| {{przyklad|Wyznaczanie miejsc zerowych wielomianu||
| |
| | |
| Jak wiadomo, wartości własne to miejsca zerowe wielomianu charakterystycznego
| |
| macierzy <math>\displaystyle P(\lambda) = \det(A - \lambda I)</math>. Zachodzi także fakt odwrotny, to
| |
| znaczy miejsca zerowe wielomianu są wartościami pewnej macierzy, np. miejsca
| |
| zerowe wielomianu
| |
| | |
| <center><math>\displaystyle p(\lambda) = p_1 \lambda^N + \ldots + p_N \lambda + p_{N+1}
| |
| </math></center>
| |
| | |
| są wartościami własnymi m.in. macierzy stowarzyszonej,
| |
| | |
| <center><math>\displaystyle A = \beginpmatrix
| |
| -p_2/p_1 & -p_3/p_1 & \cdots & -p_{N+1}/p_1\\
| |
| 1 & & & \\
| |
| & 1 & & \\
| |
| & & \ddots & \\
| |
| & & & 1
| |
| \endpmatrix
| |
| </math></center>
| |
| | |
| Funkcja Octave'a !compan(p)! wyznacza macierz stowarzyszoną dla zadanego
| |
| wielomianu o współczynnikach w wektorze <math>\displaystyle p = [p_1,\ldots,p_N, p_{N+1}]^T</math>. Z tej
| |
| macierzy korzysta następnie funkcja Octave'a !roots!, która właśnie w taki
| |
| sposób wyznacza pierwiastki wielomianów: jako wartości własne macierzy
| |
| stowarzyszonej.
| |
| }}
| |
| | |
| {{przyklad|||
| |
| Praktyczne zadanie z macierzą symetryczną
| |
| }}
| |
| | |
| W praktyce obliczeniowej spotyka się zazwyczaj kilka typów zagadnień:
| |
| * Wyznaczenie dominującej wartości własnej (to znaczy: największej co do
| |
| modułu) i odpowiadającego jej wektora własnego (a może kilku wektorów?)
| |
| * Wyznaczenie najmniejszej co do modułu wartości własnej i wektorów jej
| |
| odpowiadających (zauważmy, że to jest np. zadanie wyznaczenia ''jądra
| |
| macierzy osobliwej'' --- wtedy wiemy a priori, że szukana najmniejsza co do modułu
| |
| wartość własna to zero)
| |
| * Wyznaczenie wartości własnej najbliższej zadanej liczbie (to jest właśnie
| |
| odpowiedź na pytanie jak blisko częstości wymuszającej są częstości drgań
| |
| własnych budynku)
| |
| * Wyznaczenie wszystkich wartości własnych (na przykład, w celu znalezienia
| |
| wszystkich pierwiastków zadanego wielomianu)
| |
| * Wyznaczenie wszystkich wartości i wektorów własnych (tzw. pełne
| |
| zagadnienie własne)
| |
| | |
| Jak domyślamy się, dla macierzy rozrzedzonych dużego wymiaru pełne zagadnienie
| |
| własne jest zbyt kosztowne, gdyż najczęściej macierz wektorów własnych --- nawet
| |
| dla macierzy rzadkiej --- jest gęsta.
| |
|
| |
|
| Ponieważ w zastosowaniach bardzo często pojawiają się macierze rzeczywiste
| | <center><math>x_k\,=\,\phi( x_{k-1}),\qquad k\ge 1</math></center> |
| symetryczne (powyższe przykłady pokazują, że nie tylko!) szczegółową analizę
| |
| metod numerycznych ograniczymy do tego przypadku, gdyż wtedy zachodzi
| |
|
| |
|
| {{twierdzenie|o symetrycznym zadaniu włanym|| | | {{twierdzenie|Banacha, o zbieżności iteracji prostej|| |
|
| |
|
| Każda macierz rzeczywista symetryczna <math>\displaystyle A</math> wymiaru <math>\displaystyle N</math> ma rozkład
| | Niech <math>D_0</math> będzie domkniętym |
| | podzbiorem dziedziny <math>D</math>, |
|
| |
|
| <center><math>\displaystyle A = Q\Lambda Q^T, | | <center><math>\overline D_0\,=\,D_0\subset D, |
| </math></center> | | </math></center> |
|
| |
|
| gdzie <math>\displaystyle Q\in R^{N\times N}</math> jest ortogonalna (tzn. <math>\displaystyle Q^TQ = I</math>), a jej kolumnami są
| | w którym <math>\phi</math> jest odwzorowaniem zwężającym. |
| wektory własne <math>\displaystyle A</math>, natomiast <math>\displaystyle \Lambda\in
| | To znaczy, <math>\phi(D_0)\subset D_0</math>, oraz istnieje stała |
| R^N</math> jest diagonalna z
| | <math>0\le L<1</math> taka, że |
| wartościami własnymi <math>\displaystyle A</math> na diagonali:
| |
|
| |
|
| <center><math>\displaystyle \Lambda = \beginpmatrix \lambda_1 & & \\ & \ddots & \\ & & | | <center><math>\|\phi( x)-\phi( y)\|\,\le\,L\,\| x- y\|, |
| \lambda_N\endpmatrix . | | \qquad\forall x, y\in D_0. |
| </math></center> | | </math></center> |
|
| |
|
| }}
| | Wtedy równanie |
| | |
| ===Uwarunkowanie zadania===
| |
|
| |
|
| {{twierdzenie|Bauer-Fike||
| | <center><math> |
| | | x\,=\,\phi( x). |
| Niech <math>\displaystyle A\in R^{N\times N}</math> będzie diagonalizowalna, to
| |
| znaczy dla pewnej macierzy <math>\displaystyle X</math> zachodzi
| |
| | |
| <center><math>\displaystyle X^{-1} A X = \beginpmatrix \lambda_1 & & \\ & \ddots & \\ & &
| |
| \lambda_N\endpmatrix ,
| |
| </math></center> | | </math></center> |
|
| |
|
| a więc (gdyż macierz po prawej stronie jest podobna do <math>\displaystyle A</math>) <math>\displaystyle \lambda_i\in C</math>,
| | ma dokładnie jedno |
| <math>\displaystyle i=1,\ldots,N</math> są
| | rozwiązanie <math>x^*</math>, oraz |
| wartościami własnymi <math>\displaystyle A</math>. Rozważmy macierz zaburzoną <math>\displaystyle \tilde{A}</math> i jakąś jej
| |
| wartość własną <math>\displaystyle \tilde{\lambda}</math>. Wtedy istnieje wartość własna <math>\displaystyle \lambda_j</math>
| |
| macierzy <math>\displaystyle A</math> taka, że
| |
|
| |
|
| <center><math>\displaystyle |\lambda_j - \tilde{\lambda}| \leq </math> cond <math>\displaystyle _2(X) ||A - \tilde{A}||_2. | | <center><math>x^*\,=\,\lim_{k\to\infty} x_k, |
| </math></center> | | </math></center> |
|
| |
|
| | dla dowolnego przybliżenia początkowego |
| | <math>x_0\in D_0</math>. |
| }} | | }} |
|
| |
|
| Ponieważ dla rzeczywistej macierzy symetrycznej macierz przejścia <math>\displaystyle X</math> jest
| | {{dowod||| |
| ortogonalna,
| | Wobec |
| <math>\displaystyle X^{-1} = X^T</math>, to mamy cond <math>\displaystyle _2(X) = 1</math> i w konsekwencji zachodzi
| |
|
| |
|
| {{wniosek|Wartości własne macierzy symetrycznej są doskonale uwarunkowane|| | | <center><math>\begin{matrix} \| x_k- x_{k-1}\| & = & |
| | \|\phi( x_{k-1})-\phi( x_{k-2})\| |
| | & \le & L\| x_{k-1}- x_{k-2}\| \\ |
| | \ &\le& \cdots\;&\le&\;L^{k-1}\| x_1- x_0\|, |
| | \end{matrix}</math></center> |
|
| |
|
| Przy oznaczeniach jak {thm:Bauer-Fike}{twierdzeniu Bauera-Fike'a}, jeśli
| | dla <math>k\ge s</math> mamy |
| dodatkowo założymy, że macierz <math>\displaystyle A</math> jest rzeczywista i symetryczna, to
| |
|
| |
|
| <center><math>\displaystyle \min_{j=1,\ldots,N}|\lambda_j - \tilde{\lambda}| \leq ||A - \tilde{A}||_2. | | <center><math>\begin{matrix} \| x_k- x_s\| |
| </math></center> | | & \le & \sum_{j=s+1}^k\| x_j- x_{j-1}\| |
| | & \le & \sum_{j=s+1}^k L^{j-1}\| x_1- x_0\| |
| | & = & L^s(1+L+\cdots+L^{k-s-1})\| x_1- x_0\| \\ \ |
| | & \le & \frac{L^s}{1-L}\| x_1- x_0\|. |
| | \end{matrix}</math></center> |
|
| |
|
| }} | | Ciąg <math>\{ x_k\}_k</math> jest więc ciągiem Cauchy'ego. |
| | Stąd istnieje granica |
| | <math>\vec\alpha=\lim_{k\to\infty} x_k</math>, która należy do |
| | <math>D_0</math>, wobec domkniętości tego zbioru. Ponieważ |
| | lipschitzowskość <math>\phi</math> implikuje jej ciągłość, |
| | mamy też |
|
| |
|
| Z drugiej strony, dla macierzy niediagonalizowalnych, uwarunkowanie wartości
| | <center><math>\phi(\vec\alpha)\,=\,\phi\Big(\lim_{k\to\infty} x_k\Big) |
| własnych może być
| | \,=\,\lim_{k\to\infty}\phi( x_k) |
| dowolnie duże, co ilustruje poniższy
| | \,=\,\lim_{k\to\infty} x_k\,=\,\vec\alpha</math>,</center> |
|
| |
|
| {{przyklad|||
| | tzn. <math>\vec\alpha</math> jest punktem stałym odwzorowania <math>\phi</math>. |
| | Dla jednoznaczności zauważmy, że jeśliby istniał |
| | drugi, różny od <math>\vec\alpha</math>, punkt stały <math>\vec\beta</math>, |
| | to mielibyśmy |
|
| |
|
| <center><math>\displaystyle A_\epsilon = \beginpmatrix a & 1 \\ \epsilon & a \endpmatrix | | <center><math>\|\vec\alpha-\vec\beta\|\,=\, |
| | \|\phi(\vec\alpha)-\phi(\vec\beta)\| |
| | \,\le\,L\,\|\vec\alpha-\vec\beta\|. |
| </math></center> | | </math></center> |
|
| |
|
| Weźmy dla uproszczenia <math>\displaystyle a=0</math>.
| | Stąd <math>1<L</math>, co jest sprzeczne z założeniem, że |
| Wartości własne <math>\displaystyle A_\epsilon</math> to zera wielomianu <math>\displaystyle p_\epsilon(\lambda) = \lambda^2 - \epsilon</math>,
| | <math>\phi</math> jest zwężająca. }} |
| zatem <math>\displaystyle \lambda_\epsilon = \pm \sqrt{\epsilon}</math> i w konsekwencji
| |
|
| |
|
| <center><math>\displaystyle |\lambda_\epsilon - \lambda_0| / ||A_\epsilon - A_0|| = \sqrt{\epsilon}/\epsilon
| | Z powyższych rozważań otrzymujemy natychmiastowy |
| \rightarrow \infty,
| | wniosek dotyczący zbieżności iteracji prostych. |
| </math></center>
| |
|
| |
|
| gdy <math>\displaystyle \epsilon \rightarrow 0^+</math>, a więc uwarunkowanie takiego zadania jest
| | {{wniosek||| |
| nieskończone: dowolnie mała zmiana macierzy powoduje zaburzenie wartości
| | Przy założeniach [[twit|Uzupe�nij: twierdzenia Banacha]], |
| własnych niewspółmiernie wielkie wobec zaburzenia danych. Dodatkowo, wartości własne i wektory własne macierzy <math>\displaystyle A</math> dla
| | metoda iteracji prostych jest zbieżna co |
| ujemnego parametru <math>\displaystyle \epsilon</math> są zespolone!
| | najmniej liniowo z ilorazem <math>L</math>, tzn. |
|
| |
|
| {eigencond.png}{Zachowanie się wartości własnych macierzy <math>\displaystyle A</math> (z
| | <center><math>\| x_k- x^*\|\,\le\,L^k\,\| x_0- x^*\|</math></center> |
| parametrem <math>\displaystyle a=1</math>) w otoczeniu <math>\displaystyle \delta = 0</math>}
| |
|
| |
|
| }} | | }} |
|
| |
|
| Bardziej spektakularny przykład pochodzi od Wilkinsona:
| | {{przyklad||| |
| | | Dla ilustracji, rozpatrzmy natępujące proste |
| {{przyklad|Perfidny wielomian Wilkinsona|| | | równanie skalarne: |
|
| |
|
| Niech
| | <center><math> |
| | | x\,=\,\cos(x), \qquad \mbox{dla} \qquad x\in D= R. |
| <center><math>\displaystyle p(\lambda) = (\lambda -1)(\lambda - 2) \cdots (\lambda - 20). | |
| </math></center> | | </math></center> |
|
| |
|
| Zmiana współczynnika przy <math>\displaystyle \lambda^{19}</math> o <math>\displaystyle 10^{-7}</math> skutkuje presunięciem niektórych
| | W tym przypadku <math>\phi(x)=\cos(x)</math>. Zauważamy, że w |
| miejsc zerowych nawet o kilka jednostek na płaszczyźnie zespolonej! Poniżej
| | przedziale <math>[0,1]</math> funkcja <math>\phi</math> jest zwężająca ze |
| pokazujemy to na numerycznym przykładzie, gdzie prócz w/w zaburzenia mamy
| | stałą |
| dodatkowo z zaburzeniami powstałymi wskutek wyznaczenia współczynników
| |
| wielomianu w arytmetyce zmiennoprzecinkowej.
| |
|
| |
|
| {wilkinson.png}{Zera oryginalnego i lekko zaburzonego perfidnego wielomianu | | <center><math>L\,=\,\max_{0\le x\le 1}|\cos'(x)|\,=\,\sin(1)\,<\,1</math></center> |
| Wilkinsona.}
| |
| | |
| Jak widzimy, zera bardzo mało zaburzonego wielomianu mogą stać się wyraźnie nie-rzeczywiste!
| |
|
| |
|
| | Stąd istnieje dokładnie jedno rozwiązanie naszego równania |
| | w przedziale <math>[0,1]</math>. Rozwiązanie to może |
| | być aproksymowane z dowolnie małym błędem przy pomocy |
| | iteracji prostych, startując z dowolnego przybliżenia |
| | początkowego <math>x_0\in [0,1]</math>. |
| }} | | }} |
|
| |
|
| Jeśli chodzi o wektory własne, ich wrażliwość na zaburzenia macierzy jest
| | Zaletą iteracji prostych jest fakt, że zbieżność |
| bardziej skomplikowana i zależy m.in. od uwarunkowania wartości własnych (czego
| | nie zależy od wymiaru <math>n</math> zadania, ale tylko od stałej |
| łatwo się domyślić) oraz od tego, jak blisko siebie leżą wartości własne.
| | Lipschitza <math>L</math> (jednak w praktyce czasem sama stała Lipschitza może zależeć od |
| | wymiaru zadania...). Metoda Banacha ma szczególne zastosowanie w |
| | przypadku, gdy funkcja <math>\phi</math> jest zwężająca na całym |
| | zbiorze <math>D</math>, tzn. <math>D_0=D</math>. Jeśli ponadto <math>D</math> ma |
| | skończoną średnicę <math>\mbox{diam} (D)</math>, to dla |
| | osiągnięcia <math>\epsilon</math>-aproksymacji zera funkcji <math>f</math> |
| | wystarczy wykonać |
|
| |
|
| ===Lokalizacja wartości własnych===
| | <center><math>k\,=\,k(\epsilon)\,=\,\Big\lceil\frac |
| | | {\log(\| x_0- x^*\|/\epsilon)}{\log(1/L)}\Big\rceil |
| Jak okaże się za chwilę, czasem warto mieć ogólne rozeznanie o tym, gdzie ''z
| | \,=\,\Big\lceil\frac |
| grubsza'' leżą wartości własne danej macierzy <math>\displaystyle A</math>. W tym celu mogą być nam
| | {\log( \mbox{diam} (D)/\epsilon)}{\log(1/L)}\Big\rceil |
| pomocne dwa fakty:
| |
| | |
| {{fakt||| | |
| Dowolna wartość własna <math>\displaystyle \lambda\in C</math> macierzy <math>\displaystyle A</math> spełnia
| |
| | |
| <center><math>\displaystyle |\lambda| \leq ||A||,
| |
| </math></center> | | </math></center> |
|
| |
|
| gdzie <math>\displaystyle ||A||</math> jest dowolną normą macierzową indukowaną przez normę wektorową.
| | iteracji, niezależnie od <math>x_0</math>. Metody zbieżne dla |
| }}
| | dowolnego przybliżenia początkowego, nazywamy |
| | ''zbieżnymi globalnie''. Obie przedstawione dotychczas metody: bisekcji i |
| | Banacha, przy rozsądnych |
| | założeniach, są zbieżne globalnie. |
|
| |
|
| Rzeczywiście, skoro istnieje wektor <math>\displaystyle x\neq 0</math> taki, że <math>\displaystyle Ax = \lambda x</math>, to stąd
| | Okazuje się, że metoda iteracji prostej może być --- w bardzo szczególnych |
| <math>\displaystyle ||Ax||/||x|| = |\lambda|</math>, więc fakt powyższy wynika już z definicji normy
| | przypadkach --- zbieżna szybciej niż liniowo. Z taką sytuacją będziemy mieli, |
| macierzy:
| | gdy korzystać będziemy z metody Newtona. |
|
| |
|
| <center><math>\displaystyle ||A|| = \max_{y\neq 0}\frac{||Ay||}{||y||} \geq ||Ax||/||x||.
| | ==Wyznaczanie wektorów i wartości własnych== |
| </math></center>
| |
|
| |
|
| Drugie twierdzenie jest równie proste w dowodzie, ale daje trochę więcej
| | {{przyklad|Moj przykład||Tekst}} |
| informacji o lokalizacji widma.
| | {{algorytm|Nie robiący nic||Leż}} |
| | |
| {{twierdzenie|Gerszgorina||
| |
| | |
| Wartości własne macierzy <math>\displaystyle A</math> leżą w sumie (teoriomnogościowej) dysków <math>\displaystyle K_i</math> na
| |
| płaszczyźnie zespolonej,
| |
| | |
| <center><math>\displaystyle K_i = \{z \in C: |z - a_{ii}| \leq \sum_{j\neq i} |a_{ij}| \}, \qquad i =
| |
| 1,\ldots N.
| |
| </math></center>
| |
| | |
| }}
| |
| | |
| {{przyklad|Koła Gerszgorina|| | |
| | |
| Niech
| |
| | |
| <center><math>\displaystyle A = \beginpmatrix
| |
| 1.08930 & 1.38209 & -1.00037 & 0.69355 & 2.32178 \\
| |
| 0.14211 & 1.74696 & 1.68440 & 0.30664 & 1.26718 \\
| |
| -0.74620 & 2.02686 & -0.68293 & 0.19684 & 0.35854 \\
| |
| 0.83517 & 0.74987 & 1.71331 & 1.09765 & -0.44321 \\
| |
| 1.02132 & -2.62155 & 0.79247 & 1.11408 & 0.48076 \\
| |
| \endpmatrix
| |
| </math></center>
| |
| | |
| {gershgorindisks.png}{Lokalizacja wartości własnych macierzy <math>\displaystyle A</math> kołami Gerszgorina oraz zgrubna
| |
| lokalizacja wewnątrz okręgu
| |
| o promieniu równym <math>\displaystyle ||A||_1</math>. Dokładne wartości własne zaznaczone trójkącikami.}
| |
| | |
| }} | |
| | |
| {{przyklad|Widmo macierzy jednowymiarowego Laplasjanu|| | |
| | |
| Norma daje:
| |
| | |
| Tw. Gerszgorina daje:
| |
| | |
| W rzeczywistości,
| |
| | |
| }}
| |
| | |
| ===Metoda potęgowa, odwrotna potęgowa, RQI===
| |
| | |
| ====Metoda potęgowa====
| |
| | |
| Przypuśćmy, że wartości własne macierzy <math>\displaystyle A\in R^{N\times N}</math> spełniają
| |
| | |
| <center><math>\displaystyle |\lambda_1| > |\lambda_2| \geq \ldots \geq |\lambda_N|,
| |
| </math></center>
| |
| | |
| (to znaczy, istnieje dokładnie jedna ''dominująca'' wartość własna macierzy
| |
| <math>\displaystyle A</math>.
| |
| | |
| Załóżmy także, że istnieje baza złożona z wektorów własnych <math>\displaystyle q_1,\ldots,q_N</math> tej
| |
| macierzy (tak jest np. dla macierzy symetrycznej na mocy
| |
| {thm:symetric-eig}{twierdzenia o własnościach symetrycznego zadania
| |
| własnego}).
| |
| | |
| Kierunek własny <math>\displaystyle q_k</math> jakiejś macierzy <math>\displaystyle A</math> ma taką własność, że poddany działaniu przekształcenia
| |
| <math>\displaystyle A</math> wydłuża się <math>\displaystyle \lambda_k</math> razy, wobec tego, dowolny wektor <math>\displaystyle x\in R^N</math> poddany
| |
| działaniu <math>\displaystyle A</math> najbardziej wydłuży się w kierunku <math>\displaystyle q_1</math>. Iterując tę procedurę,
| |
| powinniśmy dostawać w wyniku wektory, w których coraz bardziej dominuje kierunek
| |
| <math>\displaystyle q_1</math>. Formalnie, niech
| |
| | |
| <center><math>\displaystyle x = \alpha_1q_1 + \ldots + \alpha_Nq_N,
| |
| </math></center>
| |
| | |
| wtedy
| |
| | |
| <center><math>\displaystyle Ax = A \left( \sum_i \alpha_iq_i \right) = \sum_i \alpha_i A q_i
| |
| = \sum_i \alpha_i \lambda_i q_i
| |
| </math></center>
| |
|
| |
|
| i w konsekwencji
| | [[Image:MNkosztexpab.png|thumb|Wersja A]] |
|
| |
|
| <center><math>\displaystyle A^kx = \sum_i \alpha_i \lambda_i^k q_i = \lambda_1^k\left(\alpha_1q_1 +
| | [[Image:MNkosztexpab.png|thumb|50px|Wersja B]] |
| \alpha_2\left(\frac{\lambda_2}{\lambda_1}\right)^kq_2 + \ldots +
| |
| \alpha_N\left(\frac{\lambda_N}{\lambda_1}\right)^kq_N \right).
| |
| </math></center>
| |
|
| |
|
| Ponieważ z założenia, że istnieje dokładnie jedna dominująca wartość własna,
| | [[Image:MNkosztexpab.png|frame|Wersja A]] |
| <math>\displaystyle \left|\frac{\lambda_N}{\lambda_1}\right| < 1</math>, to wyrażenie w nawiasie dąży do
| |
| <math>\displaystyle \alpha_1q_1</math> i w konsekwencji wektory <math>\displaystyle x_k = A^kx</math> dążą, gdy
| |
| <math>\displaystyle k\rightarrow\infty</math>, do kierunku wektora własnego <math>\displaystyle q_1</math>, to znaczy wektora
| |
| odpowiadającego dominującej wartości własnej <math>\displaystyle A</math> (o ile tylko <math>\displaystyle \alpha_1
| |
| \neq 0</math>).
| |
|
| |
|
| Szybkość zbieżności metody potęgowej jest liniowa, o współczynniku zależnym od
| | [[Image:MNkosztexpab.png|frame|200px|Wersja B]] |
| stosunku <math>\displaystyle \lambda_2/\lambda_1|</math>. W patologicznym przypadku, gdy <math>\displaystyle |\lambda_1|
| |
| \approx |\lambda_2|</math>, może więc okazać się, że metoda praktycznie nie jest
| |
| zbieżna.
| |
|
| |
|
| W praktyce nie wyznaczamy wzorem <math>\displaystyle x_k = (A^k)\cdot x</math>, lecz raczej korzystamy z
| | <div style="margin-top:1em; padding-top,padding-bottom:1em;"> |
| metody iteracyjnej
| | <span style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie: Ciag dalszy</span> |
| | <div class="exercise"> |
|
| |
|
| [title<nowiki>=</nowiki>Metoda potęgowa]
| | Spróbuj obniżyć koszt wyznaczania <math>\exp(x)</math> dla dużych <math>x</math>! |
| <math>\displaystyle x_0</math> <nowiki>=</nowiki> dowolny wektor startowy; k <nowiki>=</nowiki> 0; | |
| while( !stop )
| |
| {
| |
| <math>\displaystyle y_k</math> <nowiki>=</nowiki> <math>\displaystyle Ax_{k-1}</math>;
| |
| <math>\displaystyle x_k</math> <nowiki>=</nowiki> <math>\displaystyle y_k/||y_k||_\infty</math>;
| |
| k++;
| |
| }
| |
|
| |
|
| Warunek normowania ma m.in. na celu zapobieżenie powstawania nadmiaru i
| | <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none"> |
| niedomiaru (gdy <math>\displaystyle |\lambda_1| < 1</math>, to <math>\displaystyle ||A^kx|| \rightarrow 0</math>, a gdy
| | <div style="font-size:smaller; background-color:#efe"> Rzecz w tym, że dla dużych <math>x</math>, trzeba wziąć bardzo dużo wyrazów w szeregu |
| <math>\displaystyle |\lambda_1| > 1</math>, to <math>\displaystyle ||A^kx|| \rightarrow \infty</math>). Przy okazji, | | Taylora. Czy można tak wykombinować, by w rezultacie wziąć ich mniej? </div> |
| <math>\displaystyle ||y_k||_\infty \rightarrow |\lambda_1|</math>, a więc mamy także sposób na wyznaczenie | | </div></div> |
| przybliżenia dominującej wartości własnej.
| |
|
| |
|
| Zazwyczaj jako warunek stopu wybiera się kryterium małej poprawki, <math>\displaystyle ||x_k -
| | </div></div> |
| x_{k-1}|| \leq \epsilon</math>, lub warunek małego residuum, <math>\displaystyle ||Ax_k - \lambda_{1,k}
| |
| x_k||\leq \epsilon</math>, gdzie <math>\displaystyle \lambda_{1,k}</math> jest przybliżeniem <math>\displaystyle \lambda_1</math>
| |
| dostępnym na <math>\displaystyle k</math>-tej iteracji.
| |
| | |
| {}{Zasada działania metody potęgowej}
| |
| | |
| Metoda potęgowa doskonale sprawdza się, gdy macierz <math>\displaystyle A</math> jest macierzą
| |
| rozrzedzoną --- np. w przypadku macierzy Google'a.
| |
| | |
| ====Odwrotna metoda potęgowa====
| |
| | |
| Zauważmy, że dla dowolnej macierzy kwadratowej <math>\displaystyle A</math> o wartościach własnych
| |
| <math>\displaystyle \lambda_k</math> i odpowiadających im wektorach własnych <math>\displaystyle q_k</math>, mamy:
| |
| * Macierz <math>\displaystyle A-\sigma I</math> ma wartości własne <math>\displaystyle \lambda_k - \sigma</math> oraz wektory
| |
| własne <math>\displaystyle q_k</math>,
| |
| * Jeśli dodatkowo <math>\displaystyle A</math> jest nieosobliwa, to macierz <math>\displaystyle A^{-1}</math> ma wartości
| |
| własne <math>\displaystyle 1/\lambda_k</math> oraz wektory własne <math>\displaystyle q_k</math>
| |
| | |
| Łącząc te dwie własności mamy, że
| |
| | |
| {{stwierdzenie|Transformacja widma macierzy||
| |
| | |
| Macierz <math>\displaystyle (A-\sigma I)^{-1}</math> (o ile istnieje),
| |
| to ma wartości własne równe <math>\displaystyle \frac{1}{\lambda_k - \sigma}</math> i wektory własne
| |
| identyczne z <math>\displaystyle A</math>.
| |
| }}
| |
| | |
| Skoro tak, to jeśli najbliższą <math>\displaystyle \sigma</math> wartością własną <math>\displaystyle A</math> jest <math>\displaystyle \lambda_j</math>,
| |
| wówczas metoda potęgowa zastosowana do macierzy <math>\displaystyle (A-\sigma I)^{-1}</math> zbiegnie do
| |
| <math>\displaystyle q_j</math>. To prowadzi do następującego algorytmu, odwrotnej metody potęgowej:
| |
| | |
| [title<nowiki>=</nowiki>Odwrotna metoda potęgowa]
| |
| <math>\displaystyle x_0</math> <nowiki>=</nowiki> dowolny wektor startowy; k <nowiki>=</nowiki> 0;
| |
| while( !stop )
| |
| {
| |
| Rozwiąż układ równań <math>\displaystyle (A-\sigma I)y_k = x_{k-1}</math>;
| |
| <math>\displaystyle x_k</math> <nowiki>=</nowiki> <math>\displaystyle y_k/||y_k||_\infty</math>;
| |
| k++;
| |
| }
| |
| | |
| ====Metoda Rayleigh====
| |
| | |
| Z własności metody potęgowej, metoda odwrotna potęgowa jest zbieżna tym
| |
| szybciej, im bliżej <math>\displaystyle \lambda_j</math> jest przesunięcie <math>\displaystyle \sigma</math> (w stosunku do
| |
| pozostałych wartości własnych). Dlatego dobrze byłoby --- dla zwiększenia
| |
| szybkości zbieżności iteracji --- poprawiać wartość przesunięcia <math>\displaystyle \sigma</math>,
| |
| korzystając z dotychczas wyznaczonego wektora <math>\displaystyle x_k \approx q_j</math> i ilorazu
| |
| Rayleigh:
| |
| | |
| <center><math>\displaystyle \lambda_j = \frac{q_j^TAq_j}{q_j^Tq_j} \approx \frac{x_k^TAx_k}{x_k^Tx_k}
| |
| </math></center>
| |
| | |
| [title<nowiki>=</nowiki>Metoda RQI (Rayleigh Quotient Iteration)]
| |
| <math>\displaystyle x_0</math> <nowiki>=</nowiki> dowolny wektor startowy; <math>\displaystyle \sigma_0</math> <nowiki>=</nowiki> przybliżenie <math>\displaystyle \lambda_j</math>; k <nowiki>=</nowiki> 0;
| |
| while( !stop )
| |
| {
| |
| Rozwiąż układ równań <math>\displaystyle (A-\sigma_k I)y_k = x_{k-1}</math>;
| |
| <math>\displaystyle x_k</math> <nowiki>=</nowiki> <math>\displaystyle y_k/||y_k||_2</math>;
| |
| <math>\displaystyle \sigma_{k+1}</math> <nowiki>=</nowiki> <math>\displaystyle x_k^TAx_k</math>;
| |
| k++;
| |
| }
| |
| | |
| (wybierając normowanie wektora <math>\displaystyle x</math> w normie euklidesowej upraszczamy co nieco
| |
| algorytm).
| |
| | |
| Wielką zaletą metody RQI jest jej szybkość zbiezności: kwadratowa gdy wartość
| |
| własna jest pojedyncza, a nawet sześcienna w przypadku macierzy symetrycznej.
| |
| | |
| Wadą metody RQI jest to, że na każdym jej kroku należy rozwiązywać układ równań
| |
| z ''inną'' macierzą.
| |
| | |
| {{uwaga|Gdy złe uwarunkowanie pomaga...||
| |
| | |
| Przez pewien czas numerycy odnosili się do tej metody z rezerwą,
| |
| twierdząc, i słusznie, że im lepszym przybliżeniem <math>\displaystyle q_j</math> będzie <math>\displaystyle \sigma_k</math>, tym
| |
| bardziej rośnie uwarunkowanie <math>\displaystyle A-\sigma_k I</math>, a tym samym --- błąd numerycznego
| |
| rozwiązywania układu z tą macierzą będzie coraz większy i metoda będzie tracić
| |
| stabilność. Tymczasem okazuje się, że --- choć rzeczywiście tak jest ---
| |
| wektor błędu ma kierunek praktycznie zgodny z kierunkiem poszukiwanego wektora
| |
| <math>\displaystyle q_j</math>, a tym samym tylko ''pomaga'' w zbieżności metody!
| |
| }}
| |
|
| |
|
| ===Metoda dziel i rządź, metoda QR=== | | <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em"> |
|
| |
|
| {}{Secular equation}
| | Bez zmieniejszenia ogólności, załóżmy, że <math>x\geq 0</math>. |
|
| |
|
| ===Biblioteki=== | | Jeśli <math>x = k + t</math>, gdzie <math>t \in [0,1)</math>, a <math>k</math> jest całkowite, to oczywiście |
| | <center><math> |
| | e^x = e^{k+t} = e^k e^t = e^k \exp(t)</math></center> |
|
| |
|
| LAPACK zawiera w sobie kolekcję doskonałych narzędzi do rozwiązywania zadania
| | Tak więc zadanie redukuje się do wyznaczenia <math>\exp(t)</math> dla ''małego'' <math>t</math> oraz |
| własnego. {WYmienic}
| | do co najwyżej <math>k</math> dodatkowych mnożeń potrzebnych do wyznaczenia całkowitej |
| | potęgi <math>e^k</math> (ile mnożeń ''naprawdę'' wystarczy?). Pamiętaj, przyjęliśmy, że |
| | znamy reprezentację numeryczną liczby <math>e</math>. |
|
| |
|
| ARPACK rozwiązuje zadanie własne dla macierzy rozrzedzonych, zwłaszcza
| | <div class="code" style="background-color:#e8e8e8; padding:1em"><pre> |
| symetrycznych.
| | [Wersja B] |
| | function [y, N] <nowiki>=</nowiki> expb(x, epsilon) |
| | k <nowiki>=</nowiki> floor(x); t <nowiki>=</nowiki> x - k; |
| | [y, N] <nowiki>=</nowiki> expa(t, epsilon); |
| | for i <nowiki>=</nowiki> 1:k |
| | y *<nowiki>=</nowiki> e; |
| | end |
| | N +<nowiki>=</nowiki> (k+2); |
| | end |
| | </pre></div> |
| | |
| | [[Image:MNkosztexpab.png|thumb|300px|Wersja B]] |
|
| |
|
| Funkcja !eig! w Octave wyznacza wszystkie wartości własne (lub pary własne)
| | </div></div></div> |
| zadaniej gęstej macierzy. Jak dotąd, tylko MATLAB potrafi skorzystać z ARPACKa
| |
| dla wyznaczenia fragmentów widma macierzy rzadkiej, za pomocą funkcji
| |
| !eigs!.
| |