MN13LAB

Zagadnienie własne

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

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

Ć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 $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

$B=\left(\begin{array}{cc}1& \frac{\sqrt{{ϵ}_{\text{mach}}}}{2}\\ \frac{\sqrt{{ϵ}_{\text{mach}}}}{2}& 1\end{array}\right)$

lub perfidny wielomian Wilkinsona. Dla macierzy $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 $Ax=\lambda x$, wybierając taką współrzędną $k$, dla której $||x|{|}_{\mathrm{\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 $O(N^2)$, oszacować odległość wartości własnych macierzy od zera. Jeśli więc w macierzy $A$ zachodzi

$|a_{ii}|>\sum _{j\ne i}|a_{ij}|,\mathrm{\forall }i,$

to oczywiście $\lambda =0$ nie może być wartością własną $A$, a skoro tak, to $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 $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 $A$ i odpowiadający jej wektor własny przy użyciu odwrotnej metody potęgowej?

A jeśli macierz $A$ jest (numerycznie) osobliwa?

Rozwiązanie

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

Jeśli macierz $A$ jest (numerycznie) osobliwa, należy zastosować odwrotną metodę potęgową z przesunięciem $\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 $O(N^3)$ do $O(N^2)$.

Ćwiczenie

Zbadaj, jak bardzo zmiana zera na $ϵ\approx 0$ w macierzy

$\left(\begin{array}{cc}a& 1\\ 0& b\end{array}\right)$

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

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

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

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \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 $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 $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 $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 $LD{L}^T$ macierzy $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ą $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ć $\sigma =\lambda (1+\sqrt{{ϵ}_{\text{mach}}})$. Startujemy z losowego wektora $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ć!