PKow: Różnice pomiędzy wersjami

Z Studia Informatyczne
Przejdź do nawigacjiPrzejdź do wyszukiwania
Przemo (dyskusja | edycje)
m Zastępowanie tekstu – „,↵</math>” na „</math>,”
 
(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ż}}
&#0025;
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!.

Aktualna wersja na dzień 21:48, 11 wrz 2023

Flash

<flash>file=Test19.swf|width=750|height=530</flash> http://osilek.mimuw.edu.pl/images/2/29/Test19.swf

Metoda iteracji prostej Banacha

&#0025; Zupełnie inne, i jak się okaże --- przy odrobinie sprytu bardzo skuteczne --- podejście do wyznaczania miejsca zerowego jest oparte na metodzie Banacha.

Najpierw nasze równanie nieliniowe

f(x)=0

przekształcamy (dobierając odpowiednią funkcję ϕ) do równania równoważnego (tzn. mającego te same rozwiązania)

x=ϕ(x).

Następnie, startując z pewnego przybliżenia początkowego x0, konstruujemy ciąg kolejnych przybliżeń xk według wzoru

xk=ϕ(xk1),k1

Twierdzenie Banacha, o zbieżności iteracji prostej

Niech D0 będzie domkniętym podzbiorem dziedziny D,

D0=D0D,

w którym ϕ jest odwzorowaniem zwężającym. To znaczy, ϕ(D0)D0, oraz istnieje stała 0L<1 taka, że

ϕ(x)ϕ(y)Lxy,x,yD0.

Wtedy równanie

x=ϕ(x).

ma dokładnie jedno rozwiązanie x*, oraz

x*=limkxk,

dla dowolnego przybliżenia początkowego x0D0.

Dowód

Wobec

xkxk1=ϕ(xk1)ϕ(xk2)Lxk1xk2 Lk1x1x0,

dla ks mamy

Parser nie mógł rozpoznać (nieznana funkcja „\begin{matrix}”): {\displaystyle \begin{matrix} \| x_k- x_s\| & \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}}

Ciąg {xk}k jest więc ciągiem Cauchy'ego. Stąd istnieje granica α=limkxk, która należy do D0, wobec domkniętości tego zbioru. Ponieważ lipschitzowskość ϕ implikuje jej ciągłość, mamy też

ϕ(α)=ϕ(limkxk)=limkϕ(xk)=limkxk=α,

tzn. α jest punktem stałym odwzorowania ϕ. Dla jednoznaczności zauważmy, że jeśliby istniał drugi, różny od α, punkt stały β, to mielibyśmy

αβ=ϕ(α)ϕ(β)Lαβ.

Stąd 1<L, co jest sprzeczne z założeniem, że

ϕ jest zwężająca.

Z powyższych rozważań otrzymujemy natychmiastowy wniosek dotyczący zbieżności iteracji prostych.

Wniosek

Przy założeniach Uzupe�nij: twierdzenia Banacha, metoda iteracji prostych jest zbieżna co najmniej liniowo z ilorazem L, tzn.

xkx*Lkx0x*

Przykład

Dla ilustracji, rozpatrzmy natępujące proste równanie skalarne:

x=cos(x),dlaxD=R.

W tym przypadku ϕ(x)=cos(x). Zauważamy, że w przedziale [0,1] funkcja ϕ jest zwężająca ze stałą

L=max0x1|cos(x)|=sin(1)<1

Stąd istnieje dokładnie jedno rozwiązanie naszego równania w przedziale [0,1]. 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 x0[0,1].

Zaletą iteracji prostych jest fakt, że zbieżność nie zależy od wymiaru n zadania, ale tylko od stałej Lipschitza L (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 ϕ jest zwężająca na całym zbiorze D, tzn. D0=D. Jeśli ponadto D ma skończoną średnicę diam(D), to dla osiągnięcia ϵ-aproksymacji zera funkcji f wystarczy wykonać

k=k(ϵ)=log(x0x*/ϵ)log(1/L)=log(diam(D)/ϵ)log(1/L)

iteracji, niezależnie od x0. 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.

Okazuje się, że metoda iteracji prostej może być --- w bardzo szczególnych przypadkach --- zbieżna szybciej niż liniowo. Z taką sytuacją będziemy mieli, gdy korzystać będziemy z metody Newtona.

Wyznaczanie wektorów i wartości własnych

Przykład Moj przykład

Tekst

Algorytm Nie robiący nic


 Leż
Wersja A
Wersja B
Wersja A
Wersja B

Ćwiczenie: Ciag dalszy

Spróbuj obniżyć koszt wyznaczania exp(x) dla dużych x!

Wskazówka
Rozwiązanie