MN13LAB

Ćwiczenia: zagadnienie własne

Ć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 Dodaj link: twierdzenie Gerszgorina.

Wskazówka

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

Ć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: Jakość kryteriów stopu w metodzie potęgowej

Porównaj ze sobą kryteria małej poprawki i małego residuum dla metody potęgowej. Jak zmieni się odpowiedź, gdy wiadomo, że dominującą wartością własną jest 1?

Rozwiązanie

gdyż wtedy

$||x_k-x_{k-1}||=||A{x}_{k-1}/||A{x}_{k-1}|{|}_{\mathrm{\infty }}-x_{k-1}||\approx 1/|{\lambda }_1|||A{x}_{k-1}-|{\lambda }_1|x_{k-1}||$

czyli wielkość poprawki mierzy wielkość residuum.

Ć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(ϵ)$.

Ć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

Dla rzeczywistej macierzy symetrycznej udowodnij następujący fakt dotyczący uwarunkowania zadania wyznaczenia wektorów własnych tej macierzy:

Ć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

Wyznacz numerycznie zera Dodaj link: wielomianu Wilkinsona korzystając z jakiejś metody rozwiązywania pełnego zagadnienia własnego i zbadaj, jak drobne zmiany w jego współczynnikach (w bazie naturalnej) wpływają na wynik.

Ć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. metodą Newtona. Zapisz wzory takiej iteracjii 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.

Ć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. potęgowej

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ć sobie, że istnieją rzeczywiste wartości własne, bierzemy macierz nieparzystego wymiaru $N$.

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 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ć!