Metody numeryczne fizyki/Wyznaczanie wektorów własnych i wartości własnych dla dowolnej macierzy

Z Wikibooks, biblioteki wolnych podręczników.
Przejdź do nawigacji Przejdź do wyszukiwania
Metody numeryczne fizyki
Metody numeryczne fizyki
Wyznaczanie wektorów własnych i wartości własnych dla dowolnej macierzy

Licencja
Autor: Mirosław Makowiecki
Absolwent UMCS Fizyki Komputerowej Uniwersytetu Marii Curie-Skłodowskiej w Lublinie
Email: miroslaw(kropka)makowiecki(małpa)gmail(kropka)pl
Dotyczy: książki, do której należy ta strona, oraz w niej zawartych stron i w nich podstron, a także w nich kolumn, wraz z zawartościami.
Użytkownika książki, do której należy ta strona, oraz w niej zawartych stron i w nich podstron, a także w nich kolumn, wraz z zawartościami nie zwalnia z odpowiedzialności prawnoautorskiej nieprzeczytanie warunków licencjonowania.
Umowa prawna: Creative Commons: uznanie autorstwa oraz miejsca pochodzenia książki i jej jakikolwiek części, a także treści, teksty, tabele, wykresy, rysunki, wzory i inne elementy oraz ich części zawarte w książce, i tą książkę, nawet w postaci przerobionej nie można umieszczać w jakikolwiek formie na czasopismach naukowych, archiwach prac, itp.
Autor tej książki dołożył wszelką staranność, aby informacje zawarte w książce były poprawne i najwyższej jakości, jednakże nie udzielana jest żadna gwarancja, czy też rękojma. Autor nie jest odpowiedzialny za wykorzystanie informacji zawarte w książce nawet jeśli wywołaby jakąś szkodę, straty w zyskach, zastoju w prowadzeniu firmy, przedsiębiorstwa lub spółki bądź utraty informacji niezależnie, czy autor (a nawet Wikibooks) został powiadomiony o możliwości wystąpienie szkód. Informacje zawarte w książce mogą być wykorzystane tylko na własną odpowiedzialność.


Poznamy tutaj metody obliczania wektorów i wartości własnych dla macierzy o n wierszach i n kolumnach, czyli kwadratowej, a także poznamy metody obliczania wartości własnych dla tej macierzy znając definicję algebry, który można poznać ze standardowego kursu algebry.

Wstęp do obliczania wartości i wektorów własnych dla macierzy kwadratowej[edytuj]

  • Wartości i wektory własne dla macierzy kwadratowej o stopniu n, czyli , oznaczamy kolejno przez x, λ, wtedy te wielkości spełniają warunek Axx.
  • Wartości własne liczymy ze wzoru det(AI)=0, który przedstawia wielomian charakterystyczny, który przyrównamy tu do zera.
  • Gdy rozpatrzymy macierz transponowaną AT, to jego wektory i wartości własne nazywamy lewostronnymi wartościami i wektorami własnymi. Widmem macierzy , czyli Sp(A) nazywamy zbiór wartości własnych λ12,...,λn. Widma macierzy kwadratowej jego odpowiednika transponowanego są takie same.
  • Jeśli przez λL oznaczymy wartości własne lewostronne, a przez λP wartości własne macierzy kwadratowej A, to wektory własne macierzy A i macierz transponowanej AT są takie, że są do siebie prostopadłe według definicji iloczynu skalarnego, a dowód jego przebiega:<xL,Axp>=λp(xL,xP), a także zachodzi <xL,Axp>=<ATxL,xPL<xL,xP>, stąd odejmując te dwa ostatnie równania otrzymujemy (λLp)<xL,xP>=0.
  • Jeśli wartości własne macierzy A są λ1, λ2,..,λn, to wartościami własnymi dla wielomianu p(t) przy t=A są liczby p(λ1), p(λ2),...,p(λn), to z definicji wspomnianego wielomianu wektory własne macierzy p(A) są takie same jak wektory własne macierzy A.
  • Gdy λ jest wielomianem charakterystycznym macierzy A, który przyrównamy do zera, to A musi być macierzą zerową.
  • Gdy zdefiniujemy macierz P, który jest macierzą podobieństwa i zachodzi P-1AP=B, to macierze A i B są do siebie podobne.
  • Macierz i m≥n jest macierzą ortogonalną, gdy zachodzi warunek .
  • Kolumny i wiersze macierzy ortogonalnej są wektorami do siebie ortogonalnymi.
  • Macierz diagonalna D jest macierzą podobną do macierzy A, gdy istnieje macierz ortogonalna Q, która jest macierzą ortogonalnie podobną, taką że zachodzi QTAQ=D.
  • Macierze symetryczne A mają wektory własne o składowych rzeczywistych i rzeczywiste wartości własne.
  • Macierz diagonalna D jest podobna do macierzy A na w sposób D=P-1AP, które są do siebie podobne i mają takie same wartości własne zespolone.
  • Każda macierz trójkątna, która jest w ogólności macierzą zespoloną, jest macierzą unitarnie podobną do dowolnej macierzy kwadratowej.

Zdefiniujmy macierz zwaną macierzą Jordana rzędu k, w której λ jest w ogólności liczbą zespoloną taką, że:

(6.1)
  • Jeśli macierze J1,J2,...,Jr będziemy nazywać klatkami Jordana (6.1), wtedy każdą macierz A przy definicji macierzy podobnej P, która w ogólności jest macierzą zespoloną, możemy przedstawić w postaci:
(6.2)
  • Kwadrat normy euklidesowej dla normy ||⋅||E jest większy lub równy sumie kwadratów modułów wartości własnych macierzy A.
(6.3)
  • Widmo macierzy, które jest sumą macierzy A i skalaru pomnożonej przez macierz jednostkową jest równa sumie widma wspomnianej macierzy i skalaru c.
(6.4)
której powyższe wyrażenie jest zbiorem, którego elementy są sumą liczb λi i stałej "c", czyli λ1+c,λ2+c,...,λn+c.
  • Koło domknięte o środku aii i promieniu równej sumie elementów na wierszu bez elementów na diagonali spełnia warunek:
(6.5)
to wtedy widmo macierzy A zawiera się sumie zbiorów (6.5) dla wszystkich i od i=1 do n. Jeśli stworzymy zbiór dla k rozłącznych kół określonych według powyższego schematu, to w tak utworzonym zbiorze leży k wartości własnych należącej do macierzy A.
  • Nierówność ρ(A)≥ρ(B) zachodzi, gdy moduły elementów macierzy B są mniejsze lub równe elementom macierzy A.
  • Jeśli macierz A ma wartości własne λ, to macierz B=A-cx><x ma dokładnie takie same wartości i wektory własne jak macierz A dla wektorów własnych różnych od x, a dla wektora własnego równego x odpowiada jemu wartość własna równa λ-c.

Błędy zaokrągleń dowolnej macierzy wraz z jego wartościami i wektorami własnymi[edytuj]

Będziemy się tutaj zajmować błędami zaokrągleń macierzy A i jego wartości własnej, gdy błędy zaokrągleń naszej macierzy dążą do zera, czyli dla błędy zaokrągleń wartości własnej dążą do zera, czyli dla macierze dokładnej przy zanikającym błędach zaokrągleń. Podajmy kilka definicji dotyczącej odległości pomiędzy zbiorami i odległości pomiędzy liczbą a zbiorem.

  • Odległości pomiędzy widmem macierzy A a liczbą "z", czyli zbiorem wartości własnych Sp(A)={λ12,...,λn} nazywamy liczbę zdefiniowanej w postaci , które przedstawimy według oznaczenia dist(z-Sp(A)).
  • Dla dwóch widm macierzy A i macierzy B, które są macierzami kwadratowymi, to odległość pomiędzy nimi przy zdefiniowanej wartości własnej dla tychże dwóch wspomnianych macierzy λi i γip jest dla kolejnych wskaźników i=1,2,..,n.
  • Jeśli błędy zaokrągleń macierzy A dążą do zera, to widmo Sp(AA) dąży do widma Sp(A), wtedy powiemy ||Sp(A)-Sp(A)||→0.
  • Macierz A i macierz diagonalna D, jeśli są podobne do siebie, to zachodzi warunek P-1AP=D, wtedy odległość wartości własnej macierzy przybliżonej A i macierzy dokładnej AA jest napisana poprzez oszacowanie dist(λ,Sp(A))≤α, przy której definicję stałej α jest opisywana poprzez
    α=||P-1||2||P||2||δA||2.
  • Wskaźnik uwarunkowania macierzy A, która służy według zadania do wyliczania wartości własnej macierzy A podobnej do diagonalnej nazywamy liczbę K2=||P-1||2||P||2.
  • Przy macierzy przekształcenia P i dowolnej macierzy symetrycznej A jest zadaniem dobrze uwarunkowanym, gdy zachodzi K2(P)=1.
  • Przekształcenie jednej macierzy w drugą nazywamy operację matematyczną w postaci formuły B=RAR-1, która to macierz A jest macierzą nieosobliwa, co może to pogorszyć lub polepszyć uwarunkowanie wspomnianej macierzy, bo zachodzi K2(P)=||P-1||2||P||2, a dla macierzy drugiej podobnej do poprzedniej uwarunkowanie macierzy jest równe K2=||(RP)-1||2||RP||2.
  • Gdy mamy przekształcenie B=QAQT przy macierzy ortogonalnej Q, to ono nie zmienia uwarunkowana wartości własnych macierzy, tzn. są spełnione warunki ||QP||2=||P2|| i ||(QP)-1||2=||P-1Q-1||2=||P-1||2.
  • Weźmy sobie dwie macierze A i D, która ta ostatnia macierz jest macierzą diagonalną, a wektory xi i yi dla i=1,...,n są odpowiednio wektorami własnymi i wektorami własnymi lewostronnymi wspomnianej macierzy, i oznaczając wartość własną macierzy AA przez γi, wtedy możemy napisać:
(6.6)

wtedy poszczególne składniki w sumie (6.6) przepisujemy według dwóch schematów:

(6.7)
(6.8)

Dla dowolnego zaburzenia dla dowolnej macierzy widmo macierzy A γ∈Sp(AA) przy przekształceniu P-1AP=J jest tak napisane, by było spełnione oszacowanie przy definicji α=K2(P)||δA||2, a także napiszemy wzór dla dostatecznie małych zaburzeń macierzy A poprzez oszacowanie:

(6.9)
(6.10)

Przejdźmy teraz do dowodu wzoru (6.10), które przedstawimy przy pomocy macierzy Jordana. Weźmy sobie teraz dowolną macierz A, która jest podobna do macierzy Jordana J według przekształcenia P-1AP=J. Weźmy sobie wartość własną macierzy AA, czyli γ, którego równanie własne przepisujemy wedle schematu (AA)xx, wtedy równanie własne macierzy Jordana jest ( J+P-1AP)yy, jeśli przyjmować będziemy E=P-1A)P, wtedy możemy napisać (JI)y=-Ey, wtedy lewą stronę ostatnio wspomnianego równania przepisujemy w postaci rozwiniętej:

(6.11)

Poszczególne macierze występujące w macierzy (6.11), czyli JmIm, możemy przepisać według schematu:

(6.12)

Macierz odwrotną do macierzy prostej JI piszemy wychodząc od wzoru (6.12), tzn.:

(6.13)

Macierz odwrotną do macierzy (6.11) przy definicji poszczególnych jego kratek (6.13) piszemy według:

(6.14)

Ze wzoru powyżej podanego możemy napisać równość y=-(JI)-1Ey, wtedy z definicji normy możemy napisać nierówność, która wynika z ostatniej równości 1≤||(JI)-1||2||E2, wtedy na podstawie tego wyrażenie (6.14) i piszemy:

(6.15)

Jeśli oznaczymy ρ=dist(γ,Sp(A)), wtedy wyrażenie (6.14) przepisujemy wedle:

(6.16)

Jeśli oznaczymy przez ρk=k||E||2, wtedy na pewno dla ρ<1 powiemy , a jeśli ρ>1, wtedy możemy napisać 1+ρ+..+ρk≤kρk-1, co za tym idzie ρ<kα, wtedy na podstawie tego otrzymujemy ostateczny wzór (6.9). Gdy ρ jest dostatecznie małe, tzn. zachodzi warunek ρ<1, powiemy:

(6.17)

Możemy połączyć wzory (6.17) ze wzorem (6.15) i w ostateczności mamy (6.10). Weźmy sobie wartości własne λ12,...,λm dla macierzy A, a także mamy wartości własne γi dla macierzy dokładnej AA. Gdy zaburzenie dąży do zera, to wartości własne macierzy dokładnej i zaokrąglonej dążą do siebie.

  • Niech wektor własny odpowiadający macierzy zaokrąglonej oznaczymy przez xi, a wektor własny macierzy dokładnej oznaczmy przez yi, wtedy minimalna odległość pomiędzy tymi wektorami oznaczymy przez:
(6.18)

Jeśli oznaczymy zbiór wartości własnych przez dim Xj=1, wtedy iloczyn skalarny pomiędzy nimi określamy poprzez kąt θ pomiędzy tymi wektorami, który możemy wyliczyć ze wzoru <yi,xi>/||xi||⋅||yi||=sinθ.

Twierdzenie
Weźmy sobie macierz symetryczną y, a także jej wartości własne λi, wtedy odległość pomiędzy wartościami własnymi określamy jako:, wtedy dla zaokrągleń wspomnianej macierzy δA, dla którego zachodzi ||δA||2 dla której to macierzy dokładnej AA odpowiadają wartości własne γi, wtedy odległość (6.10) pomiędzy wektorami własnymi macierzy zaokrąglonej i dokładnej jest równa:
(6.19)
Dowód
Równanie własne macierzy dokładnej możemy przepisać wedle schematu (AA)yiiyi, co po przekształceniu mamy: (Aii)yi=-δAyi. Możemy napisać dolne oszacowanie pomiędzy wartościami własnymi macierzy dokładnej i zaokrąglonej |γii|≥||δA||2, wtedy na podstawie tego piszemy ||(AiI)y||2≥(d-||δA||2)||||yii||, co na podstawie tego powiemy: ||yi||2(d-||δA||2)≥||(AiI)y⊥||2≥||yi||2 ≥||A||2). Jeśli oznaczymy przez ρ(yi,xi), co na podstawie tego możemy napisać (6.19).
  • Weźmy sobie macierz symetryczną A, która jest symetryczna, a także możemy napisać ||r||=ρ, to wtedy znajduje się co najmniej jedna wartość własna w przedziale . Dla wartości własnej macierzy, które znajdują się poza tym przedziałem o odległośc co najmniej d>ρ, wtedy rysujemy , co jest błędem wektora własnego obliczonego zgodnie z (6.18), przy którym Xj=X jest podprzestrzenią rozpiętą na wektorach własnych macierzy A, które odpowiadają wartością własnym z wyżej podanego przedziału.

Korzystając z definicji równania własnego, w którym poszczególne obliczenia są dokonane dokładnie, czyli dla wyrażenia Axx=0, wtedy wektor reszt obliczamy, gdy poszczególne obliczenia w ostatnio wspomnianym równaniu są dokonane w przybliżeniem (zaokrągleniem) sposobem:


(6.20)

W obliczeniach (6.20) dokonaliśmy przybliżenia . Będziemy przyjmować, że długość słowa binarnego jest bardzo duża, tzn. zachodzi (n+2)⋅2-t<0,1, wtedy na podstawie tego napiszemy:

(6.21)

Jeśli będziemy przyjmować , i jeśli będziemy przyjmować, że , wtedy na podstawie rozważań przedstawionym w tym zdaniu i wzoru (6.21) możemy napisać:

(6.22)

Gdy posiadamy zaokraglone wartości i i obliczenia dokonujemy z użyciem podwójnej precyzji, wtedy dowiadujemy się przy oszacowaniu błędu przy liczeniu wyrażenia Axx, że:

(6.23)

Znajdowanie miejsc zerowych dowolnej macierzy[edytuj]

W przypadku dowolnej macierzy wartości własne leżą w kole Ci o środku w punkcie aii i mającej promień równy promieniowi spektralnemu ρ(A). Dla dowolnej macierzy A zachodzi ρ(A)≤||A||p przy którym p=1,2,∞,E. Wartości własne zatem należą do zbioru {z:|z|≤||A||p}. Gdy p=∞ to oszacowanie wartości własnych piszemy przez |z|≤||A||, ale ponieważ zachodzi Sp(A)=Sp(AT), to również możemy powiedzieć |λi|≤||A||1.

Wykorzystanie metod potęgowych przy wyznaczaniu poszczególnych wartości własnych i wektorów własnych dla dowolnej macierzy[edytuj]

Metoda potęgowa polega na wyznaczaniu wartości i wektorów własnych dowolnej macierzy A. Ta metoda jest bardzo skuteczna. Algorytm ten polega na obraniu dla i=0 dowolnego wektora początkowego xi, którego norma dla ∞ jest równa jeden, tzn. ||x0|| i dalej ustalamy dalszy krok iteracji. Następnie obliczamy vi+1=Axi, a także mi+1=||vi+1||. Jesli mi+1=0 przerywamy krok iteracji, jeśli nie, to dalszym krokiem iteracji jest:

(6.24)

Jeszli dalej i+1≥ITM, to wtedy algorytm zatrzymujemy, jeśli nie, to przechodzimy do etapu pierwszego liczenia wektora vi+1. Okazuje się, że ciąg wektorów x0,x2,...,x2j jest ciągiem zbieżnym do wektora x, który jest wektorem własnym równania własnego macierzy A, również ciąg {mi} jest zbieżny do pewnej wartości m. Jeśli zachodzi ||xj+1-xj||→0, to wtedy mamy Ax=mx, a jeśli w przeciwnym wypadku ||xj+1-xj||→2||x||, to wtedy zachodzi Ax=-mx. Dla tych przypadków λ=m lub λ=-m znajdujemy wektor własny w wyniku powyższych iteracji dla wektora x, którego norma jest równa jeden dla p=∞. W ogólnym przypadku dla uzyskanego ciągu wektorów własnych x2i dla macierzy A, to może być on rozbieżny. Gdy jednak istnieje przekształcenie P, dla której macierz A jest podobna do macierzy diagonalnej λ12,..,λn, to wtedy możemy napisać twierdzenie:

Twierdzenie
Ogólnie rzecz mówiąc, jeśli mamy wartości własne λ1, λ2,..,λn, które są wartościami własnymi macierzy A, w których poszczególne elementy są różne, ale mają równe za to moduły, wtedy przy dowolnym wyborze elementu startowego x0 uzyskany ciąg wektorów x0, x2, x4,... jest ciągiem zbieżnym do pewnego wektora będącego wektorem własnym macierzy A.
Dowód
Dowód przeprowadźmy dla przypadku macierzy symetrycznej A, której wektory własne są wektorami bazy ortogonalnej w przestrzeni Rn, której wektory bazy oznaczymy przez e1, e2, e3,...,en, której odpowiadają wartości własne |λ1|≥|λ2|≥...≥|λn|. Oznaczmy wartość początkową wektora iteracji przez x01e12e2+...+αnen. Może się zdarzyć, że wszystkie współczynniki liniowości są równe zero oprócz jednej lub wartość własna może być wielokrotna. Weźmy sobie λ=0, wtedy możemy obrać tak x0, by było można napisać Ax0=0 dla ||x0||=1. Natomiast, gdy λ≠0, to wtedy przyjmując oznaczenia

aiei+...+αi+kei+k, δ0i+k+1ei+k+1+...+αnen z powodów oczywistych otrzymujemy:



(6.25)

Przyjmując oznaczenia, co do δ0, wówczas możemy napisać z założenia, że wektor δ0 dąży do zera, a także że wartości własne macierzy A nie rosną wraz z rosnącym wskaźnikiem:


(6.26)
Według obliczeń (6.25) i (6.26) od razu dochodzimy do wniosku, że zachodzą dwie tożsamości, tzn. oraz , wtedy na podstawie tego powiemy , że jest zawsze zbieżny do , co jest wektorem własnym macierzy A dla wartości własnej λ. Oznaczmy wzór na yL w postaci yL=xL(sqn λ), i wtedy mamy dla λ>0, że yj=xj, i dla λ<0 mamy y2j+1=-x2j+1. Dla pierwszego przypadku mamy ||xi+1-xi||→0, a dla drugiego przypadku dostajemy ||xi+1-xi||→2||x||.

Wykorzystanie metod Hessenberga przy wyznaczaniu poszczególnych wartości własnych dla dowolnej macierzy[edytuj]

Każdą macierz Hassenberga nazywamy macierz zapisaną w postaci sumy dwóch macierzy, tzn. macierzy trójdiagonalnej (5.93) T i macierzy trójkątnej górnej U, zdefiniowanej według H=T+U. Wiadomo jednak hij=0 dla j<i-1. Jeśli jednak potrafimy wyznaczyć wartości własne macierzy H, to również potrafimy wyznaczyć wartości własne dowolnej macierzy. Według prostego algorytmy QR każdą macierz A możemy rozłożyć na iloczyn macierzy ortogonalnej Q i macierzy trójkątnej R, którą zapisujemy według schematu A=QR. Gdy macierz A jest nieosobliwa, to poszczególne kolumnach macierzy Q możemy otrzymać dokonując ortogonalizacji macierzy A metodą Grama-Schmidta, wtedy kolumny macierzy R są zbudowane z współczynników rozwiniecie z ortogonalizowanej macierzy A. Zbudujmy teraz iterację, która jest algorytmem QR, która opisuje pewien ciąg macierzy H=H(i),H(i),..., która jest zbudowana w postaci dwóch wzorów:

(6.27)
(6.28)

We wzorach (6.27) i (6.28) wszystkie macierze Q(i), a także H dla i=1,2,... są na pewno macierzami Hessenerga. Ogólnie rzecz biorąc powyższe dwa wzory możemy połączyć w jedną całość jako:

(6.29)

Jeśli mamy wartości własne λ12,...,λn, które są parami sprzężone, tzn. mającej równe moduły, wtedy nie wszystkie elementy pod diagonalą dążą do zera. Macierz Hessenberga można otrzymać z ciągu macierzy H dla wskaźnika górnego górnego nieskończonego, wtedy otrzymujemy np. macierz:

(6.30)

We wzorze (6.30) gwiazdki (*) oznaczają elementy zbieżne, a kropki (⋅) oznaczają wartości dążące do ustalonych wartości. Mówiąc ogólnie macierz Hi nie musi mieć ogólnie postaci (6.30), ale może mieć za to postać:

(6.31)

Mając macierz (6.31) wykorzystując przy tym wzór (6.29), to w wyniku końcowych obliczeń otrzymujemy tożsamości H(1)=H2=...=H(i) przy macierzy H ortogonalnej, to w rozkładzie H=QR macierz R występująca w tymże ostatnim wzorze jest macierzą jednostkową. Zbieżność elementów występująca na diagonali w macierzy Hessenberga jest często powolna, aby temu zaradzić możemy dokonać przesunięcia wartości własnych, w ten sposób tworzymy ciąg macierzy H(1),H(2),H(3),..., które definiujemy przy pomocy dwóch wzorów:

(6.32)
(6.33)

Łącząc wzory (6.32) i (6.33) w jedną całość dostajemy końcowy iteracyjny wzór, który generuje elementy zbieżne do macierzy Hessenberga, którego schemat iteracyjny przepisujemy w postaci:

(6.34)

Gdy obierzemy odpowiednio ki to otrzymujemy szybszą zbieżność do macierzy Hessenberga H(i), współczynniki ki dobiera je się w takiej postaci by one były równe wartością własnym kolejnych elementu ciągu macierzy Hessenberga, tzn. , gdzie elementy ki i ki+1 są równe kolejnym elementom macierzy kwadratowej dwuwymiarowej usytuowanej w H(i) w prawym dolnym rogu tejże wspomnianej macierzy. Ale te elementy mogą być również zespolone, wieć rachunki prowadzą do rachunku na macierzach zespolonych, wtedy należy wykonać dwie dalsze iteracje, ale ta metoda zawodzi, gdy mamy doczynienia z macierzami np. (6.31).

Znajdowanie wartości własnych dowolnej macierzy poprzez doprowadzenie jej do postaci Hessenberga[edytuj]

Istnieje wiele metod sprowadzania dowolnej macierzy A do macierzy Hessenberga poprzez wykorzystanie wiadomości o podobieństwie dwóch macierzy przy przekształceniu P. Istnieją trzy takie metody, tzn. metoda Householdera, metoda Givensa i na samym końcu metoda eliminacji Gaussa. Pierwsze dwie metody potrzebują i mnożeń w arytmetyce zmiennoprzecinkowej, dzięki której otrzymujemy macierz Hessenberga H, której błąd zaokrąglenia jest ||E||E≤k2εn2||A||E dla pierwszej metody oraz ||E||E≤k1εn1/2||A||E dla drugiej metody. Stałe k1 i k2 występujące w dwóch ostatnich wzorach są oczywiście rzędu jedynki, które zależą od sposobu zaokrąglania w komputerze, czyli na maszynie cyfrowej. Natomiast liczba mnożeń w metodzie w eliminacji Gaussa jest , którego ilość operacji jest istotnie mniejsza niż w poprzednich metodach, i ta metoda jest powszechnie stosowana w obliczeniach numerycznych przy wyznaczaniu wartości własnych dowolnej macierzy A. W tej ostatniej metodzie wykonujemy n-2 przekształceń uzyskując macierze A=A(1),A(2),...,A(n-1)=H, gdzie ta ostatnia macierz jest macierzą Hessenberga. Macierz A(i) w której i-1 początkowych kolumn są również i-1 kolumnami macierzy Hessenberga H. Przekształcenie macierzy A ze wskaźnika górnego i do wskaźnika górnego i+1 postępujemy w taki sposób, by wybierać element o największym module z elementów , a gdy tychże elementów jest więcej, to wybieramy element pierwszy z nich, a jeśli mamy doczynienia z zerowymi elementami, to powinno być A(i)=A(i+1), dalej wyznaczamy element A(i+2). Jeśli element o największym module znajduje się na k-tym wierszu dla k≥i, to należy przedstawić k-ty wiersz i i+1-tą kolumną. Dalej liczymy wielkość dla j=i+2,i+3,..,n, które zgodnie z wyborem elementu maksymalnego zachodzi , co dzięki temu nasza metoda jest numerycznie stabilna. Następnie i+1-ty wiersz odejmujemy od j-tego wiersza pomnożonej przez i j-tą kolumnę dodajemy od i+1-tej kolumny pomnożonej przez dla j=1+2,i+3,..,n, wtedy przy takim postępowaniu uzyskana macierz A(i+1) jest podobna do macierzy A(i) przy elementach , które przy naszym postępowaniu są równe zero, co oznacza, że i-ta kolumna A(i+1) jest taka sama jak i-ta kolumna macierzy Hessenberga H.

Przegląd metod wyznaczania wektorów własnych metodą QR,LR i metodą iteracji odwrotnej[edytuj]

Niech macierz P jest macierzą przekształcenia, w którym macierz A jest podobna do macierzy B, jeśli potrafimy wyznaczyć wektory własne według równania Bxx, to wtedy wektor y=Bx jest wektorem własnym macierzy A. Jeśli założymy, że macierz B jest macierzą trójkątną górną, ale posiada różne wartości własne, wtedy wartość własna macierzy B odpowiada liczba λi=bii, a wektorowi własnemu wektor x(i), które te wielkości wyliczamy z formuł:

(6.35)

Wyznaczanie wektorów własnych metodą QR[edytuj]

W tej metodzie należy wyznaczyć macierze Q(1),Q(i),..., a następnie znając macierz przekształcenia P wyznaczyć macierz Hessenberga H z macierzy A według P-1AP=H. Dalej ze wzorów (6.27) i (6.28) lub ze wzorów (6.32) i (6.33) możemy wyznaczyć macierze H(i+1)=Q(i)TH(i)Q(i), a także wyznaczamy kolejne macierze pomocnicze Z(i+1)=Z(i)Q(i) z warunkiem początkowym Z(1)=P. Dalej musimy stwierdzić, czy macierz H(N) ma elementy pod diagonalą dostatecznie bliskie zero, jeśli tak, to możemy przystąpić do wyznaczania wektorów własnych macierzy Hessenberga ze wzoru:

(6.36)

Wyznaczanie wektorów własnych metodą LR[edytuj]

Ta metoda jest bardzo podobna do poprzednio opisywanej metody, tutaj należy wyznaczyć ciąg macierzy A=A(1),A(2),.., które liczymy z formuł dla i=1,2,..., które to L i U są odpowiednio macierzami trójkątnymi dolnymi i górnymi:

(6.37)
(6.38)

Kolejne macierze A(i) możemy wyznaczyć ze wzoru A(i+1)=(L(i))-1A(i)L(i)m, które są podobne do macierzy , który ten wzór wynika ze wzorów (6.37) i (6.38), gdzie macierze A(i) dążą do macierzy trójkątnej. W celu wyznaczenia wektorów własnych dla macierzy A trzeba wyznaczyć obiekty Z(i+1)=Z(i)L(i) i Z(1)=I. Aby otrzymać rozkład A(i)=L(i)U(i) należy wykonać mnożeń, czyli trzy razy mniej niż przy rozkładzie A(i)=QR(i), co sugeruje przewagę tejże metody nad metoda poprzednią, jednak stosując poprzednią metodę wykonujemy mniejszy błąd zaokrągleń, ale natomiast gdy macierz A jest macierzą symetryczną i dodatnio określoną to ta metoda jest stabilna numeryczna i szybsza od poprzednio opisywanej metody. Podobnie jak w poprzedniej metodzie stosuje się tutaj przesunięcie A(i)-kiI, mając początkowe macierze A sprowadzane się je do macierzy Hessenberga, co zmniejsza liczbę obliczeń przy rozkładzie macierzy Hessenberga H=LU przy nakładzie pracy .

Wyznaczanie wektorów własnych metodą iteracji odwrotnej[edytuj]

Metody QR, LR są bardzo podobne do metody potęgowej, dlatego że w nim następuje znikanie niektórych kolumn macierzy A na podstawie jego różnych wartości własnych, ale zbieżność w metodzie potęgowej jest jednak bardzo powolna, ale za to w tychże metodach przy właściwym wyborze przesunięcia ki, która jest bliska wartości własnej λn, która jej dobrym przybliżeniem jest element następuje zwiększenie zbieżności przy wyznaczaniu wektorów własnych. Szybkość zbieżności jest zależna od , następuje wtedy zwiększenie szybkości zbieżności tejże metody. Wartościami własnymi macierzy A są zaś liczby λ1, λn, ...,λn takie, że spełniają nierówności słabe |λ1|≥|λ2|≥....≥|λn|, to stąd wynika, że wartościami własnymi macierzy odwrotnej do A-kiI są liczby , wtedy można tak dobrać ki, by uzyskać zbieżność wektorów zależnych od wskaźnika górnego "i" do wektora własnego macierzy A, czyli x z omówionym powyżej współczynnikiem. Omówmy metodę iteracji potęgowej. Weźmy i=0 i obierzmy taki wektor x0 by był spełniony warunek ||x0||=1 i ustalamy maksymalną wartość iteracji ITM. Dalej wyznaczamy wektor vi+1=(A-kiI)-1xi i skalar mi+1=||vi+1||, a także również wektor xi+1=vi+1/mi+1. Jeśli i+1≥ITM zatrzymujemy algorytm, a jeśli nie, to wtedy poprzedni krok wykonujemy dla i:=i+1. Liczby m1, m2,... są stosowane do wyznaczania wartości własnej λn, tak jak w metodzie potęgowej. Wartość ki jest przybliżeniem wartości λn, gdy jednak jest ona dostatecznie bliska do niej, wtedy macierz (A-kiI) jest oczywiście macierzą źle uwarunkowaną, co prowadzi do dużych błędów numerycznych przy wyznaczaniu wektora vi+1, a jeśli nasz wektor własny jest dobrze uwarunkowany, wtedy zaokrąglenia przy wyznaczaniu układu równań zapisanej w postaci macierzowej (A-kii)vi+1=xi nie będą zbyt duże, bo musimy się zabezpieczyć się by wektor vi+1 był taki by jego norma była dostatecznie duża. Weźmy sobie kin+en, by potem można było wyznaczyć wektor vi+1, w której macierz A jest zapisana z pewnym zaokrągleniem, który zapisujemy w postaci układu równań zapisanej w postaci macierzowej (A+E-kiI)vi+1=xi przy założeniu ||E||=e. Napiszmy:



(6.39)

Przy opisie (6.39) należy stosować taki układ równań w postaci macierzowej (A-kiI)vi+1=xi, by norma błędu ||E|| była mała, ale duży współczynnik uwarunkowania (A-kiI) nie jest wtedy bardzo groźny, stąd wynika, że wektor reszt (An)xi+1 jest bardzo mały dla dobrze uwarunkowanego wektora x, stąd wynika, że dobrym przybliżeniem jest do niego xi+1.

Ogólne metody rozkładu dowolnej macierzy na iloczyn macierzy QR[edytuj]

Podamy tutaj twierdzenia dotyczące rozkładu macierzy A na iloczyn macierzy QR, a później omówimy metodę rozkładu dowolnej macierzy Grama-Schmidta, Householdera i Givensa.

Twierdzenie
Dowolną macierz przy założeniu m≥n na k≤n niezależnych w nim liniowo kolumn, wtedy mamy rozkład A=QR, przy czym macierz jest macierzą ortogonalną, a macierz jest macierzą trapezoidalną górną, tzn. posiadająca elementy zerowe pod diagonalą (rij=0 dla i≤j), a jeśli mamy doczynienia z macierzą kwadratową A, to macierz R jest macierzą trójkątną górną. Należy pamiętać, że dowolną macierz możemy rozłożyć w sposób jednoznaczny na iloczyn A=QR, to wtedy elementy na diagonali są zawsze dodatnie.

Metoda rozkładu dowolnej macierzy A na iloczyn macierzy Q i R metodą Grama-Schmidta[edytuj]

Poznamy tu rozkład dowolnej macierzy A na rozkład macierzy QR metodą ortogonalizacji Gramma-Schmidta, które zastosujemy dla naszej macierzy dla kolumn a1, a2,...,n. Obierzmy sobie macierze A(1), A(2),...,A(n), które posiadają 1,2,...,.n kolumn macierzy A. Jeśli dodatkowo założymy, że pierwsza kolumna naszej macierzy jest niezerowa, to wtedy napiszmy wzory na Q(1) i na R(1), które wykorzystamy do liczenia macierzy A(1):

(6.40)
(6.41)
(6.42)

Jeśli mamy już macierze początkowe, tzn. (6.40), (6.41) i (6.42), to wtedy możemy policzyć następne elementy tychże macierzy, tzn. Q(i+1), Ri+1, a na jej podstawie policzyć macierz A(i+1):

(6.43)
(6.44)
(6.45)

We wzorach (6.43), (6.44) i (6.45) występują pewne nieznane wielkości, które to wyznaczymy wedle poniższych formuł dla i=1,2,...,n-1:

(6.46)
(6.47)
(6.48)

Jeśli si+1=0, to zachodzi równość Q(i+1)=Q(i), a jeśli mamy si+1≠0, to wtedy możemy napisać qi+1=pi+1/si+1. Łatwo udowodnić, że macierz (6.44) jest macierzą trapezoidalną, a macierz Q(i) jest macierzą ortogonalną, a macierz A(i) spełnia zależność (6.45). Wyniku użycia metody Grama-Schmidta jest możliwa utrata ortogonalności macierzy Q(i) dla i=1,2,..,n w wyniku obliczeń numerycznych, które powstają wyniku zaokrąglenia wyników obliczeń. W tym celu należy przeprowadzić reortogonalizację pi+1, ale to zwiększa nakład obliczeń numerycznych. Przy rozkładzie macierzy A=QR, która z założenia jest macierzą nieososbliwą należy wykonać dodawań i i "n" pierwiastkowań, które wykonujemy bez ponownej reorogonalizacji.

Metoda rozkładu dowolnej macierzy A na iloczyn macierzy Q i R metodą Householdera[edytuj]

Weźmy sobie macierz kwadratową H o k wierszach i k kolumnach przy zdefiniowanym wektorze u=z-α||z||2e1 i definicji liczby τ=||z||22-α(z,e1):

(6.49)

Macierz zdefiniowana w punkcie (6.49) nazywamy macierzą Householdera. Podziałajmy macierzą (6.49) na wektor z=[z1,z2,...,zk]T przy wyżej podanych definicjach u i τ:

(6.50)
Aby pominąc błędy zaokrągleń spowodowane dzieleniem, która podowują nadmierny wzrost błędów numerycznych należy tak obrać odpowiednią liczbę τ, wtedy wprowadza się definicję stałej α w postaci: α=-sqn z1, co powoduje α=-1 dla z1≥0 i α=1 dla z1<0. Macierz przekształcenia P(1)=H(1) ustala pierwszą kolumnę macierzy A zapisaną w postaci schematu: do macierzy przekształcenia. Zdefiniujmy sobie macierz powstałej z iloczynu lewostronnego macierzy przekształceń o wskaźnikach i=1,...,n-1 i macierzy A wedle definicji P(n-1)P(n-2)....P(2)P(1)A=R. Zdefiniujmy macierz przekształcenia o wskaźniku górnym dwa jako:, wtedy macierz H(2) stworzona jest na podstawie wierszy i kolumn macierzy P(1)A o numerach i=2,..,n, w której pierwszą kolumną jest , którą jeśli podziałamy macierzą H, to wtedy mamy . Mając , to możemy podstawić za z do równania (6.49) i w ten sposób otrzymać . Idąc dalej
(6.51)

wtedy macierz H(n-1) stworzona jest na podstawie wierszy i kolumn macierzy P(n-2)...P(2)A o numerach n-1 i n, w której pierwszą kolumną jest , którą podziałamy macierzą H, otrzymujemy , mając z podstawiamy go do (6.49), w ten sposób otrzymujemy H(n-1). Na podstawie powyższych rozważań możemy rozłożyć macierz A na iloczyn macierzy Q i R, czyli A=QR, wtedy poszczególne czynniki w tym iloczynie piszemy według jego definicji Q=P(1)P(2)⋅...⋅P(n-1) oraz R=QTA.

Metoda rozkładu dowolnej macierzy A na iloczyn macierzy Q i R metodą Givensa[edytuj]

Metoda Givensa jest bardzo podobna do metody Householdera. Weźmy sobie wektor , a wynik działania macierzy Givensa na wektor z przedstawiamy α||z||2⋅[1,0,...,0]T, tzn. dla α=±1 mamy Gz=α||z||2e1. Macierz Givensa definiujemy sposobem:

(6.52)

W macierzy (6.52) są tak odpowiednio dobrane liczby c(j), s(j), by podczas działania macierzy G(j) na wektor pionowy , by j+1 współrzędna wyniku działania była równa zero. Opisane tutaj liczby wyznaczamy według schematów:

(6.53)
(6.54)
(6.55)
(6.56)

Napiszmy macierz G jako iloczyn macierzy Givensa o wskaźnikach górnych i=2,...,n, którego definicja jest G=G(2)⋅...⋅G(n-1)G(n). Podczas działania macierzy G(n) na wektor pionowy n wyzerowuje jego element zn, i tak idąc dalej do końca macierz G(2), który podczas działania na wektor z wyzerowuje element z2 do zera. Utwórzmy sobie macierz , który podczas działania na wektor (jest to pierwsza kolumna macierzy A) otrzymujemy wynik . Weźmy sobie macierz G2, który to działamy na wektor pionowy (jest to pierwsza kolumna utworzonej macierzy o n-1 wierszach i kolumnach położonej w prawym dolnym rogu macierzy G1A) otrzymujemy wynik . Na samym końcu zbudujmy sobie macierz Gn-1 podczas działania na wektor , który jest macierzą o wymiarze 2x2 utworzoną z wierszy i kolumn macierzy Gn-2Gn-1...GA otrzymujemy wektor . Podczas rozkładu macierzy A na iloczyn Q i R otrzymujemy przy dokonanych obliczeniach przy powyższym schemacie otrzymujemy następujące ich definicje , a . Liczba operacji potrzebna do wyznaczenia macierzy Givensa i macierzy R jest w postaci mnożeń i dodawań i pierwiastkowań.

Wyznaczenie wartości własnej dla symetrycznej trójdiagonalnej macierzy[edytuj]

Macierz trójdiagonalna jest definiowana wedle schematu (5.93), a dodatkowo jeśli jest symetryczna, to jego definicja jest:

(6.57)

Omówimy tutaj dwie metody, tzn. metodę bisekcji i metodę QR jako metody wyznaczania odpowiednich wartości własnych.

Wyznaczanie wartości własnych dla trójdiagonalnej macierzy symetrycznej metodą bisekcji[edytuj]

Weźmy sobie dowolną liczbę λ, wtedy możemy napisać wielomian charakterystyczny ω(λ)=det(AI). Ponieważ co do założenia macierzy A wielomian ω(λ) możemy wyznaczyć ze związków rekurencyjnych zapisanej w postaci:




(6.58)
  • Jeśli wszystkie elementy a2, a3,...,an są niezerowymi wartościami, to macierz (6.57) ma wszystkie wartości własne różne od siebie, a także na podstawie (6.57) można wygenerować ciąg wartości zależnych od λ, tzn. ω0(λ), ω1(λ),...,ωn(λ), które dla ωi(λ), gdy i<n, spełniają warunki ωi-1(λ)ωi+1(λ)<0, lub gdy spełniony jest warunek ωn(λ)≠0, to liczba zmian znaków opisywanego powyżej ciągu, dla i=1,..,n, jest równa liczbie wartości własnych macierzy (6.58). A ponadto gdy ωn(λ)=0, to λ jest wartością własną macierzy (6.58), i tyle jest wartości własnych oczywiście mniejszych niż λ, ile nastąpiło zmian znaków ciągu dla naszego powyżej opisywanego ciągu dla i=0,..,n-1.

Powyższa metod chociaż jest metodą stabilną, ale poszczególne elementy ciągu ω0(λ), ω1(λ),...,ωn(λ) mogą osiągać zbyt duże wartości, nawet dla małych liczb naturalnych "n", co następuje, gdy wartość λ odbiega zbytnio od wartości własnych macierzy (6.57), z tego względu ciąg generowany przez algorytm (6.58) przekształcamy w ciąg pi: , wtedy jego poszczególne elementy są generowane przez ciąg wynikającego z przekształconego algorytmu (6.58), jako:

(6.59)
(6.60)

Gdy ciąg generowany (6.60) ma element pi-1 równy zerowy, to ten nasz algorytm musi ulec modyfikacji. Wyznaczmy dokładną wartość wyrażenia pi(λ), który to przepisujemy znając zaokrągloną wartość bi, ai i pi-1, na w sposób:

(6.61)

Liczba zmian znaków elementu dokładnego fl(pi(λ)) jest równa liczbie zmian znaków (6.61) podzielonej przez (1+ε2)(1+ε3), bo podzielenie przez to wyrażenie nie zmienia znaku powyższego wyrażenia, co wtedy tak otrzymane wyrażenie, a właściwie jej prawą stronę piszemy przez:

(6.62)

Oszacowania, co do błędów poszczególnych elementów (6.57), które powstają podczas zaokrąglania, są napisane jako:

(6.63)
(6.64)

Ze wzorów (6.63) i (6.64) wynika ||δT||<2,06ε||T||. Weźmy sobie początkowy przedział , który dzielimy pokolei na m części i obierzemy przez wartość t jako środek ostatniego przedziału i mając wartość własną λk możemy napisać jego oszacowanie:

(6.65)

Jeśli dodatkowo przyjmować będziemy β(0)=||T|| i α(0)=-||T||, to wzór (6.65) na podstawie oszacowania δT mamy w postaci:

(6.66)

Jeśli wyznaczymy wartości własne opisane metodą powyżej, to możemy przystąpić do liczenia poszczególnych składowych wektora własnego macierzy (6.57) według sposobu:

(6.67)
(6.68)
(6.69)

Uzyskiwanie rozkładu macierzy trójdiagonalnej symetrycznej metodą QR[edytuj]

Weźmy sobie macierz trójdiagonalną T (6.57), którą będziemy rozkładać na iloczyn dwóch macierzy Q i R, wtedy możemy napisać iloczyn macierzy trójdiagonalnej P(1), P(2),...,P(n-1) i macierzy A, który jest macierzą R, zapisywanej jako P(n-1)P(n-2)⋅...⋅P(2)P(1)T=R. Jeśli podziałamy macierzą P(1) na macierz trójdiagonalną T, to ona ma za zadanie wyeliminować w macierzy T element t21 do zera. Macierzą P(1) może być zarówno macierzą Householdera H(1) (6.49), jak i macierzą Givensa G(2). Wektor u jest tak zdefiniowany by miał dwie początkowe składowe. Obie te macierze, tzn. macierz Householdera lub Givensa, możemy zapisać znając elementy "a" i "b" w formie:

(6.70)

którym składowe element a i b w obu metodach odpowiednio różnią się znakiem. A więc metody Householdera i Givensa są na pewno równoważne w tym przypadku. Dla metody Givensa macierze R i Q są zapisywane jako: R=G(n)G(n-1)⋅...⋅G(3)T=R i Q=G(2)⋅...⋅G(n), która ta ostatnia macierz jest macierzą Hessenberga. Jeśli wymnożymy wszystkie czynniki przez siebie w macierzy Q i w R, to w ostateczności otrzymamy macierze w jego pełnych formach:

(6.71)
(6.72)

Poszczególne elementy pełnej macierzy Q, tzn. pi, βi i si możemy zdefiniować wedle schematu w formie:

1)
2)
3)

Odwróćmy teraz kolejność mnożenia w przedstawianiu macierzy T , wtedy otrzymujemy inną macierz , ale ona powinna być podobna ortogonalnie do macierzy T (6.57) według schematu RQ=QTTQ, wtedy na podstawie tego przedstawienie tej naszej macierzy:

(6.73)

Podczas rozkładu macierzy T (6.57) na iloczyn QR należy wykonać M=9n-9 mnożeń, D=3n-3 dodawań i n-1 pierwiastkowań. Dodatkowo, aby obliczyć macierz wedle jego schematu (6.72) należy wykonać M=4n-1 mnożeń i D=n dodawań. Ilość pamięci jakie są zużywane do wyznaczenia macierzy R I g(2),....,G(n) jest 4n-3 komórek pamięci maszyny cyfrowej.

Sprowadzanie dowolnej macierzy symetrycznych do postaci trójdiagonalnej[edytuj]

Okazuje się, że można sprowadzić dowolną macierz A do postaci trójdiagonalnej symetrycznej, tzn. istnieje taka macierz przekształcenia P, dzięki któremu jest to możliwe.

Wyznaczanie macierzy trójdiagonalnej symetrycznej metodą Householdera[edytuj]

Obierzmy sobie macierze P(1), P(2),...,P(n-2). Pierwszą macierz tego ciągu macierzy przekształceń piszemy przez:

(6.74)

Właściwości macierzy H(i) zostały omówione w rozdziale Metoda rozkładu dowolnej macierzy A na iloczyn macierzy Q i R metodą Householdera. Zdefiniujmy macierz M(2), który wychodzi z macierzy A, którą będziemy doprowadzać do postaci diagonalnej:

(6.75)

Obierzmy sobie drugą macierz przekształcenia zapisanej przy pomocy macierzy H(2) obraną w formie:

(6.76)

wtedy możemy obliczyć macierz M(3) otrzymaną z macierzy M(2) poprzez macierz przekształcenia P(2) otrzymywaną na w sposób:

(6.77)

Gdy wykonamy n-2 przekształceń macierzy A otrzymujemy w końcu M(n-1)=T, czyli macierz trójdiagonalną symetryczną. Macierze H(i) są macierzami Householdera zapisnych w postaci macierzy (6.49). Aby wyznaczyć macierz H(i)AH(i) należy policzyć wektor , a także . Wykorzystując definicję macierzy Householdera (6.49) możemy napisać:





(6.78)

Obliczenia przeprowadzone w punkcie (6.78) bardzo są potrzebne do obliczania pokolei macierzy M(i). Aby wyznaczyć trójdiagonalną symetryczną macierz należy dokonać mnożeń i dodawań oraz n-2 pierwiastkowań. Błąd w tej metodzie wyznaczania macierzy dokładnej T w jego przybliżonej wersji jest:

(6.79)

Stała k1 we wzorze (6.79) zależy od sposobu zaokrąglania na typowej maszynie cyfrowej i dla typowej maszyny cyfrowej jej wartość jest k1=20.

Wyznaczanie macierzy trójdiagonalnej symetrycznej metodą Givensa[edytuj]

Przy wyznaczaniu macierzy P(i), która jest macierzą przekształcenia, musimy skorzystać z macierzy Givensa opracowanej w punkcie Metoda rozkładu dowolnej macierzy A na iloczyn macierzy Q i R metodą Givensa, i na jej podstawie mamy macierze przekształcenia, co w rezultacie mamy ciąg M(2),M(3),...,M(n-1)=T. W tej metodzie dokonaliśmy mnożeń i dodawań, a także pierwiastkowań. Błąd w tej metodzie wyznaczania macierzy dokładnej T w jego przybliżonej wersji jest:

(6.80)