Akademia górniczo-hutnicza



Pobieranie 0,79 Mb.
Strona3/5
Data24.02.2019
Rozmiar0,79 Mb.
1   2   3   4   5

Właściwości macierzy
Macierz dodatnio określona (ang. „positive definite matrix”)

Macierz A nazywa się dodatnio określoną jeśli:



  • jest symetryczna

  • dla każdego niezerowego wektora z wartościami rzeczywistymi x

zachodzi xTAx > 0
Jest to macierz analogiczna do dodatniej liczby rzeczywistej.

Przykłady





M0 jest dodatnio określona. Dla wektora z=[z0 z1] (z0 z1 są liczbami rzeczywistymi) z definicji



co jest zawsze większe od zera.


b)

M1 nie jest dodatnio określona, ponieważ gdy z=[1 -1] iloczyn zTM1z jest równy -2.


Macierz jest również dodatnio określona jeśli wszystkie jej wartości własne są dodatnie.
Układy równań w których występuje macierz dodatnio określona mogą być rozwiązane w sposób wydajny przy pomocy metody rozkładu Cholesky’ego. Metoda ta pozwala nawet na dwukrotnie szybsze obliczenia w porównaniu do metody Gaussa, przy jednoczesnej oszczędności pamięci. Poprawiona jest także stabilność numeryczna.
Dla takiej macierzy niepotrzebny jest także wybór elementów podstawowych (jak również zmiana koleności wierszy) gdy używa się metody Gaussa-Crouta.
Inne własności:

  • macierz dodatnio określona jest zawsze odwracalna. Jej odwrotność również jest określona dodatnio.

  • suma macierzy dodatnio określonych również jest dodatnio określona


Metoda eliminacji Gaussa
Jest to jedna z najczęściej używanych metod numerycznego rozwiązywania układów równań liniowych. Może też posłużyć m.in do obliczania rzędu macierzy oraz znajdywania macierzy odwrotnej. Używa podstawowych operacji na wierszach by zredukować macierz do postaci schodkowej.
Etap 1 (etap eliminacji zmiennych)
Dany jest układ równań A(1)x=b(1), postaci:
a(1)11 x1 + a(1)12 x2 + … + a(1)1n xn = b(1)1

a(1)21 x1 + a(1)22 x2 + … + a(1)2n xn = b(1)2



………………………………………………………………………

a(1)n1 x1 + a(1)n2 x2 + … + a(1)nn xn = b(1)n


Odejmując od i-tego wiersza układu ( i=1,2,…n ) wiersz pierwszy pomnożony przez a(1)i1/a(1)11 otrzymuje się układ A(2)x=b(2) ,postaci:

a(1)11 x1 + a(1)12 x2 + … + a(1)1n xn = b(1)1

………. a(1)22 x2 + … + a(1)2n xn = b(1)2

………………………………………………………………………

………. a(1)n2 x2 + … + a(1)nn xn = b(1)n


Po tym kroku wyeliminowano niewiadomą x1 z równań leżących w wierszach o numerach i>=2. Przez analogiczne odejmowanie od wierszy układu kolejno zwielokrotnionych wierszy j>=2 otrzymuje się przekształcone układy równań A(3)x=b(3) , A(4)x=b(4), … . Po (n-1) takich eliminacjach uzyskuje się trójkątny układ równań.
Etap 2

Trójkątny układ równań z etapu 1 rozwiązuje się metodą podstawiania .


Żeby rozwiązać metodą eliminacji Gaussa układ n równań z n niewiadomymi maszyna cyfrowa musi wykonać 1/3n3 + n2 – 1/2n mnożeń i dzieleń oraz 1/3 n3 + 1/2n2 -5/6n dodawań. Złożoność algorytmu w związku z powyższym to O(n3) .
Metoda ta jest stabilna dla macierzy przekątniowo dominujących i dodatnio określonych.
Rozkład LU

Rozkład LU (ang „LU decomposition” lub „LU factorization”) to jedna z metod rozkładu macierzy kwadratowej na iloczyn. Dla macierzy A o n wierszach i n kolumnach rozkład zwraca dwie macierze kwadratowe trójkątne L oraz U o wymiarach takich samych jak wymiar A. Macierze te posiadają własność A = LU. L to macierz trójkątna dolna (tzn niezerowe elementy znajdują się jedynie na przekątnej i pod nią), zaś U to macierz trójkątna górna (sytuacja odwrotna niż w L, elementy niezerowe znajdują się na i nad przekątną).



Przykład


Zastosowanie rozkładu LU

  1. Rozwiązywanie układów równań – podstawowe zastosowanie metody w tym projekcie. Znając rozkład macierzy A możliwe jest bardzo szybkie (ze złożonością O(n2) ) znalezienie rozwiązania używając metody podstawiania (ang „forward/backward substitution”). Żeby uzyskać jednak rozkład należy użyć metody o złożoności porównywalnej do metody eliminacji Gaussa. Dlatego rozwiązywanie układów równań przy pomocy rozkładu LU jest wydajne głównie gdy potrzebne jest obliczenie wielu układów równań z tą samą macierzą A, różniących się wektorem b.

Sposób rozwiązywania układu równań wygląda następująco:



  1. Ax=b

  2. Ponieważ A=LU , układ równań przyjmuje postać LUx = b

  3. Rozwiązać równanie Ly = b ,gdzie y=Ux

  4. Rozwiązać równanie Ux = y




  1. Obliczanie macierzy odwrotnej

  2. Obliczanie wyznacznika macierzy – znając rozkład LU można bardzo szybko obliczyć wartość wyznacznika macierzy A. Wyznacznik A jest równy iloczynowi wyznaczników macierzy L i U. Ponieważ L i U są trójkątne, ich wyznaczniki równe są iloczynowi elementów na przekątnych.

det(A) = det(L)*det(U)=∏i=1..n Lii * ∏i=1..n Uii


Do przeprowadzenia rozkładu LU można użyć metody eliminacji Gaussa, lub tzw metody Doolittle’a. Obie wymagają 2n3/3 operacji zmiennoprzecinkowych. Łączna liczba działań potrzebnych do rozwiązania układu równań (łącznie z rozkładem ) jest taka sama jak w metodzie Gaussa.
Opis metody Doolittle’a

Dla i = 1,2,…,n obliczamy

uij = aij - ∑k=1..(i-1) Lik*ukj ; j = i,i+1,…,n

Lij = (aij - ∑k=1…(i-1) Ljk*uki)/uii ; j=i+1, i+2, …, n


Stabilność metody Gaussa i Doolittle’a

Metody Gaussa i Doolittle’a w przedstawionej powyżej formie nie są niezawodne. W trakcie wykonywania algorytmu występuje dzielenie przez elementy na przekątnej. Jeśli na przekątnej znajdzie się zero, proces zostanie przerwany. Rozwiązaniem tego problemu jest modyfikacja metod przez tzw. wybór elementu podstawowego (ang „pivoting”). Elementem podstawowym nazywa się element macierzy A za pomocą którego eliminuje się zmienną z dalszych równań. Elementem podstawowym w pierwotnych wariantach metody Gaussa i Doolittle’a jest właśnie element na przekątnej. Przez zmianę kolejności wierszy macierzy A umieszcza się na przekątnej element niezerowy z tej samej kolumny (częściowy wybór elementu podstawowego, ang „partial pivoting”). Algorytm wybiera element o największej wartości bezwzględnej. Ma to na celu zagwarantowanie numerycznej stabilności rozwiązania (dzielenie przez większe liczby zmniejsza błędy). Istnieją też warianty poszukujące elementu podstawowego w całej macierzy, czyli dopuszczając zamianę kolumn (całkowity wybór elementu podstawowego, ang „complete pivoting”). Stosowane są jednak rzadko ponieważ nie gwarantują poprawy stabilności, jednocześnie wprowadzając dodatkowe problemy obliczeniowe.

Zmiana elementu podstawowego nie jest konieczna dla niektórych typów macierzy.Do takich przypadków należą:


  1. macierze równocześnie nieosobliwe i diagonalnie dominujące

  2. macierze symetryczne i dodatnio określone

Rozkład LDLT

Innym rozkładem który można użyć do przyspieszenia obliczeń jest rozkład macierzy A na macierz L-macierz trókątna dolna zawierająca 1 na przekątnej, U-macierz trójkątna górna zawierająca 1 na przekątnej oraz D-macierz diagonalna o elementach na przekątnej równych elementom na przekątnej macierzy U z rozkładu LU. Można go uzyskać m.in. przez eliminację Gaussa. Rozkład taki jest jednoznaczny. Jeśli macierz jest symetryczna to U=LT, więc ostatecznie rozkład przyjmuje formę A = LDLT. Macierze L i D uzyskuje się z poniższych wzorów:


d1 =a11

Lij = (aij - ∑k=1…(j-1) cik * Ljk ) / dj ; i=2,3,…,n

cij = dj *Lij ; j=1,2,…,(i-1)

di = aii - ∑ k=1…(i-1) cik * Lik


Tą metodą możemy dla macierzy symetrycznej znaleźć rozwiązanie x przy pomocy blisko dwa razy mniejszej liczbie działań niż przy metodzie Gaussa .
Rozkład Cholesky’ego

Rozkład Cholesky’ego polega na zamianie macierzy A na iloczyn macierzy kwadratowych L i LT. Macierz A musi być symetryczna i dodatnio określona. Rozkład A = LLT nie jest jednoznaczny, chyba że liczby na przekątnej są dodatnie. Otrzymuje się wtedy następujące wzory:



; i=1,2,…,n

; j=i+1, i+2,…,n
Rozkład Cholesky’ego jest ok. 2 razy szybszy od rozkładu LU. Dodatkową zaletą jest fakt, że niepotrzebne jest stosowanie metody wyboru elementu podstawowego, ponieważ rozkład ten wymaga macierzy zdefiniowanych dodatnio.
Metody iteracyjne rozwiązywania układów równań

W metodach iteracyjnych proces rozwiązywania polega na uzyskiwaniu kolejnych przybliżeń wektora wynikowego x. Początkowa wartość x może być w niektórych metodach dowolna.

W kolejnych krokach rozwiązania uzyskuje się wyniki zmierzające do wyników rzeczywistych (pod warunkiem zbieżności metody). Metody iteracyjne mogą służyć zarówno jako jedyna metoda obliczeń, jak również jako sposób zmniejszenia błędu uzyskanego przy obliczeniach innymi algorytmami. Metody iteracyjne są często używane do rozwiązywania układów w których występuje znaczna liczba niewiadomych.

Głównymi problemami związanymi z zastosowaniem metod iteracyjnych są:



  1. konieczność uzyskania zbieżności rozwiązania –istnieją układy równań dla których kolejne kroki iteracyjne zwiększają początkowy błąd, zamiast przybliżać do poprawnego rozwiązania

  2. wybór liczby iteracji –algorytm musi zdecydować kiedy przerwać obliczenia, czyli kiedy szacowany błąd osiągnął akceptowalną przez użytkownika wartość

Wykorzystuje się wiele metod iteracyjnych. Poniżej opisane zostaną jedne z najpopularniejszych.

Metoda Jacobiego

W metodzie Jacobiego macierz A rozkłada się na składniki L+D+U, gdzie L jest macierzą poddiagonalną, D macierzą diagonalna, zaś U to macierz naddiagonalna.


A L D U
Ogólny wzór na (i+1) rozwiązanie to:

Dx(i+1) = -(L+U)x(i) + b ; i>=0

Żeby móc jednak zastosować powyższy wzór macierz wejściowa A musi mieć elementy niezerowe na przekątnej. Przed użyciem metody może być więc konieczna zamiana kolejności wierszy. Jeśli w każdej kolumnie macierzy A znajduje się taka sama liczba elementów zerowych (często spotykana sytuacja dla macierzy rzadkich),to dla uzyskania macierzy z niezerową przekątną wystarczy postępować podobnie jak przy wariantach metody Gaussa z wyborem elementu podstawowego. Nie gwarantuje to jednak zbieżności metody. Metoda nie jest także zbieżna dla każdej dodatnio określonej macierzy symetrycznej. Zbieżność jest zapewniona gdy A jest silnie diagonalnie dominująca.


Metoda Gaussa-Seidla

Metoda ta jest podobna do metody Jacobiego.Ogólny wzór to:

Dxi+1 = -Lxi+1 – Uxi +b ;i>=0

Jak widać ze wzoru, w tej metodzie uzyskuje się kolejne elementy wektora x na podstawie już obliczonych. Jest to ważne ze względu na wymagania pamięciowe, ponieważ wystarczy przechowywać jeden wektor zawierający aktualne rozwiązanie(w przeciwieństwie do metody Jacobiego). Na początku kalkulacji przyjmuje się wektor x o dowolnych wartościach (przy czym im bliższy jest on dokładnemu rozwiązaniu tym szybciej uzyskuje się rozwiązanie).

Aby móc stosować metodę Gauss-Seidla, należy przyjąć taką kolejność równań, żeby elementy na przekątnej były niezerowe. Jeśli A jest silnie diagonalnie dominująca to metoda jest zbieżna. Gdy macierz A jest symetryczna, wówczas zbieżność jest zapewniona gdy jest jednocześnie dodatnio określona(co jest kolejną zaletą nad metodę Jacobiego). Poważnym ograniczeniem stosowalności jest fakt, że obliczenia nie mogą być wykonywane równolegle(ponieważ nowe elementy wektora wynikowego powstają na podstawie już wyliczonych).
LAPACK

LAPACK (Linear Algebra PACKage, pakiet algebry liniowej) to biblioteka programistyczna do obliczeń numerycznych. Udostępnia procedury umożliwiające m.in rozwiązywanie układów równań, zagadnienia związane z wartościami własnymi oraz faktoryzację( rozkład) macierzy. Napisany został w języku Fortran, ale istnieje wiele implementacji na inne języki programowania (jak C++,Java,C# i Haskell). Do obliczeń wykorzystuje system BLAS (Basic Linear Algebra Subprograms, podprogramy do podstawowej algebry liniowej), który zapewnia wysoką wydajność na wielu platformach. Biblioteka jest dostępna na licencji BSD.

W ramach niniejszego projektu napisano prosty interfejs biblioteki LAPACK w języku C++, przystosowany do użycia obiektowego. Dokładny opis zamieszczono w dalszej części pracy.
Projekt modułu Solvera
Opis wymagań

Celem projektu było stworzenie projektu modułu Solvera do biblioteki MES. Jest wiele kryteriów które decydują o możliwościach modułu i w konsekwencji jego przydatności. Moduł powinien spełniać następujące wymagania:



Wymagania funkcjonalne

  • projekt powinien udostępniać funkcje implementujące najpopularniejsze i najbardziej wydajne metody numeryczne

  • udostępniać struktury danych przechowujące macierze i wektory z uwzględnieniem ich specjalnych właściwości pozwalających na wykorzystanie specjalnie dobranych metod obliczeniowych

Wymagania niefunkcjonalne

  • napisany w języku C++, zgodnie z istniejącymi standardami

  • napisany w zgodności z paradygmatem programowania zorientowanego obiektowo

  • stworzony projekt musi być łatwy do modyfikacji i przyszłego rozszerzenia

  • udostępnione funkcje numeryczne muszą być wydajne i zoptymalizowane

  • zapewniać prosty i przejrzysty interfejs

  • kompatybilny z projektem biblioteki MES tworzonym przez kolegów Piotra Opiekuna i Sławomira Muniaka

Diagram klas




1   2   3   4   5


©operacji.org 2017
wyślij wiadomość

    Strona główna