Personal tools

MN13LAB

From Studia Informatyczne


Zagadnienie własne

<<< Powrót do strony głównej przedmiotu Metody numeryczne

Oglądaj wskazówki i rozwiązania
Ukryj wskazówki i rozwiązania

Ćwiczenie

Dlaczego wartości własnych macierzy nie należy szukać jako miejsc zerowych wielomianu charakterystycznego?

Wskazówka

Pomyśl nie tylko o koszcie wyznaczenia wyznacznika dla dużych \displaystyle N.

Rozwiązanie

Przypuśćmy, że nawet w jakiś tajemny sposób wyznaczyliśmy już wielomian charakterystyczny naszej macierzy. Rzecz w tym, że miejsca zerowe wielomianu są bardzo wrażliwe na (nieuchronne) zaburzenie współczynników. Rozważ na przykład macierz

\displaystyle B = \begin{pmatrix}  1 & \frac{\sqrt{\epsilon_{\mbox{mach}}}}{2}\\ \frac{\sqrt{\epsilon_{\mbox{mach}}}}{2} & 1 \end{pmatrix}

lub perfidny wielomian Wilkinsona. Dla macierzy \displaystyle B polecenia Octave dają

octave:1> format short e

octave:2> a = sqrt(eps)/2 a = 7.4506e-09 octave:3> A = [1 a; a 1] A =

  1.0000e+00   7.4506e-09
  7.4506e-09   1.0000e+00

octave:4> eig(A)-1 ans =

  -7.4506e-09
   7.4506e-09

octave:5> roots([1 -2 1-a^2]) ans =

 1.0000e+00 + 8.2712e-09i
 1.0000e+00 - 8.2712e-09i

Z kolei w MATLABie dostajemy, trochę bardziej w zgodzie z intuicją,

>> format short e

>> a = sqrt(eps)/2

a =

  7.4506e-09

>> A =[1 a; a 1]

A =

  1.0000e+00   7.4506e-09
  7.4506e-09   1.0000e+00

>> eig(A)

ans =

  1.0000e-00
  1.0000e+00

>> eig(A)-1

ans =

 -7.4506e-09
  7.4506e-09

>> roots([1 -2 1-a^2])-1

ans =

    0
    0

Ćwiczenie

Udowodnij twierdzenie Gerszgorina.

Wskazówka

Rozpisz \displaystyle Ax = \lambda x, wybierając taką współrzędną \displaystyle k, dla której \displaystyle ||x||_\infty = |x_k|.

Wskaż, jak wykorzystać to twierdzenie do wykonania szybkiego testu, czy dana macierz jest nieosobliwa.

Rozwiązanie

Dowód twierdzenia znajdziesz w rozdziale 5.2

  • D. Kincaid, W. Cheney, Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa, 2006.

Twierdzenie Gerszgorina pozwala dość tanio, bo kosztem \displaystyle O(N^2), oszacować odległość wartości własnych macierzy od zera. Jeśli więc w macierzy \displaystyle A zachodzi

\displaystyle |a_{ii}| > \sum_{j \neq i} |a_{ij}|, \qquad \forall i,

to oczywiście \displaystyle \lambda = 0 nie może być wartością własną \displaystyle A, a skoro tak, to \displaystyle A musi być nieosobliwa. Naturalnie, jeśli powyższe kryterium nie jest spełnione, kwestia (nie)osobliwości macierzy pozostaje otwarta.

Zauważmy, że powyższe kryterium to nic innego, jak znana nam definicja macierzy diagonalnie dominującej: właśnie z twierdzenia Gerszgorina wynika, że macierz diagonalnie dominująca jest nieosobliwa.

Ćwiczenie

Czy warunek normowania wektora \displaystyle x_k jest konieczny, gdy metodę potęgową stosuje się do macierzy Google'a?

Rozwiązanie

Nie, bo dominująca wartość własna macierzy Google'a jest równa 1.

Ćwiczenie: Wyznaczanie najmniejszej wartości własnej

Jak wyznaczyć najmniejszą co do modułu wartość własną macierzy symetrycznej \displaystyle A i odpowiadający jej wektor własny przy użyciu odwrotnej metody potęgowej?

A jeśli macierz \displaystyle A jest (numerycznie) osobliwa?

Rozwiązanie

Zastosować odwrotną metodę potęgową z przesunięciem \displaystyle \sigma = 0. Kłopoty mogą wystąpić, gdy jest kilka takich wektorów własnych (na przykład, \displaystyle A jest macierzą obrotu).

Jeśli macierz \displaystyle A jest (numerycznie) osobliwa, należy zastosować odwrotną metodę potęgową z przesunięciem \displaystyle \sigma = O(\nu).


Ćwiczenie

Podaj sposób efektywnej implementacji metody odwrotnej potęgowej dla macierzy gęstych. Wykonaj ją korzystając z właściwych procedur LAPACKa (lub MATLABa).

Wskazówka

Rzecz tkwi w tym, by jak najmniej napracować się przy rozwiązywaniu

kolejnych układów

równań.

Wskazówka

W kolejnych iteracjach metody zmienia się tylko prawa strona układu równań!

Rozwiązanie

Należy wykonać rozkład LU (lub inny, stosownie do znanych własności macierzy) przed rozpoczęciem iteracji, a w trakcie iteracji wykorzystywać już tylko gotowe czynniki rozkładu --- redukując tym samym koszt pojedynczej iteracji z \displaystyle O(N^3) do \displaystyle O(N^2).


Ćwiczenie

Zbadaj, jak bardzo zmiana zera na \displaystyle \epsilon \approx 0 w macierzy

\displaystyle \begin{pmatrix}  a & 1 \\ 0 & b \end{pmatrix}

wpływa na zmianę jej wartości własnych.


Ćwiczenie: Rozwiązywanie zagadnienia własnego metodą Newtona

Parę własną \displaystyle (\lambda, x) można scharakteryzować jako rozwiązanie układu równań nieliniowych

\displaystyle \aligned Ax - \lambda x &= 0,\\ \frac{1}{2}x^Tx - 1 = 0, \endaligned

który można rozwiązać np. wielowymiarową metodą Newtona. Zapisz wzory takiej iteracji i porównaj tę metodę z metodą RQI.

Wskazówka

Generalnie, RQI jest lepsza. Metodę Newtona warto zmodyfikować tak, by na każdym jej kroku warunek normowania wektora \displaystyle x był spełniony, dzięki czemu zbieżność będzie szybsza.

Ćwiczenie

Napisz program w Octave, w którym sprawdzisz w warunkach kontrolowanego eksperymentu, że faktycznie odwrotnej metodzie potęgowej nie przeszkadza, że macierz \displaystyle A - \sigma I jest prawie osobliwa.

Rozwiązanie

Przykładowy program mógłby wyglądać następująco;

N = 51;
A = rand(N,N); 
[V, L] = eig(A);

Najpierw szykujemy sobie losową macierz, a potem wyznaczmy jej pary własne. Aby zagwarantować, że istnieją rzeczywiste wartości własne, bierzemy macierz nieparzystego wymiaru \displaystyle N (dlaczego?).

Potem definiujemy funkcję realizującą iterację odwrotną.

function x = finvit(A,x0,s)
N = size(A,1); x = x0;
	for i = 1:3
		x = (A-s*eye(N,N))\x;
		x = x/norm(x,2);
	end
end

Nie jesteśmy tu eleganccy i za każdym obrotem pętli dokonujemy ponownego rozkładu macierzy. Usprawiedliwiamy się tym, że program i tak wykona się szybciej niż wpiszemy tych kilka komend Octave, które zrealizują wstępny rozkład \displaystyle LDL^T macierzy \displaystyle A, a potem wykorzystają go w pętli.

Ponieważ będziemy korzystać z bardzo dobrego przybliżenia, wymusimy użycie

zafiksowanej liczby iteracji.

Dalej, losujemy (rzeczywistą) parę własną \displaystyle A:

reals = find(imag(diag(L))==0);
k = floor(rand(1,1)*size(reals,1))+1;
K = reals(k);
X = V(:,K); #wektor własny, wartość własna to L(K,K)

i lekko zaburzamy wartość własną, aby dostać \displaystyle \sigma = \lambda(1+\sqrt{\epsilon_{ \mbox{mach} }}). Startujemy z losowego wektora \displaystyle x_0 i wykonujemy 3 iteracje odwrotnej metody potęgowej.

s = L(K,K)*(1+sqrt(eps)); x = rand(N,1);

fprintf(stderr,'Szukamy %d-tej wartości, \n\t %e, z residuum %e\n',
	 K, L(K,K), norm(A*X - L(K,K)*X) );

x = finvit(A,x,s); l = x'*A*x;

Wyniki są faktycznie zachęcające:

octave:1> invit

Szukamy 51-tej wartości,

        2.884323e-01, z residuum 4.298045e-15

Po INVIT. Wartość własna

        2.884323e-01, z residuum 4.789183e-16.

Uwarunkowanie użytej macierzy: 4.708990e+10

Jak widać, oryginalny wynik uzyskany z procedury eig udało się nawet trochę poprawić!