1 Dwie metody rozwiązywania układów równań liniowych:Bezpośrednie Iteracyjne

2 Metody bezpośrednie Transformują układ równań na inny układ równań (prostszy do rozwiązania) Transformacjom podlegają równocześnie lewa i prawa strona równania:

3 Metody bezpośrednie Operacje wykorzystywane w przekształceniach:Zamiana dwóch równań miejscami Mnożenie równania przez niezerową stałą Mnożenie równania przez niezerową stałą i odjęcie wyniku od innego równania

4 Metody bezpośrednie

5 Układy Lx=c, Uy=d Rozwiązanie układu Lx=c metodą podstawiania do przodu Analogicznie rozwiązanie układu Uy=d (kod podany przy omawianiu dekompozycji LU)

6 Metoda eliminacji Gaussa-JordanaPodstawowa operacja: Eq.(j) – tzw. równanie (wiersz) piwotu

7 Przykład – faza eliminacjiDo roz

8 Przykład – faza eliminacji

9 Przykład – zapis macierzowyPostać „uzupełniona” (augmented)

10 Przykład – faza podstawiania wstecz

11 Eliminacja Gaussa-Jordana - algorytmk wierszy przekształconych k-ty wiersz – wiersz piwotu transformowany i-ty wiersz

12 Eliminacja Gaussa-Jordana - algorytmDo wyeliminowania wyraz Aik Trzeba podzielić wiersz piwotu przez λ=Aik/Akk i odjąć wynik od transformowanego wiersza:

13 Eliminacja Gaussa-Jordana - algorytm„Jądro” algorytmu: k - wiersz piwotu - musi przebiegać: k=1,2, ..., n-1 i - transformowany wiersz dla ustalonego wiersza piwotu k – musi przebiegać: i = k+1, k+2, ..., n j – kolumna w wierszu musi przebiegać: j = k, k+1, ..., n, takie przekształcenie jak na wierszu i - na wyrazie b(i) wektora reprezentującego prawą stronę układu równań

14 Eliminacja Gaussa-Jordana - algorytmKod źródłowy:

15 Eliminacja Gaussa-Jordana - optymalizacjePętla operacji na elementach wiersza i zastąpiona jedną operacją wektorową:

16 Eliminacja Gaussa-Jordana - optymalizacjeZastosowana uzupełniona reprezentacja macierzy A (A=[A b])

17 Eliminacja Gaussa-Jordana - optymalizacjeWyrazy pod główną przekątną nie są nam potrzebne (Lx=c), więc nie jest konieczne aby je zerować:

18 Eliminacja Gaussa-Jordana - optymalizacje

19 Eliminacja Gaussa-Jordana - optymalizacjeNie trzeba przekształcać wiersza, jeśli na odpowiedniej pozycji już jest zero:

20 Eliminacja Gaussa

21 Faza podstawiania wstecz – kodRozwiązanie ostatniego równania:

22 Faza podstawiania wstecz – kodOgólna sytuacja: jesteśmy na wysokości równania k, wszystkie równania poniżej (k+1...n) zostały już rozwiązane, mamy więc: xk+1 do xn są już znane, a więc:

23 Faza podstawiania wstecz – kod

24 Algorytm Gaussa-Jordana - kod

25 Dekompozycja LU Potrzeba: układy równań rozwiązywane wielokrotnie, dla różnych wektorów danych bi (RHS – prawe strony równań) Uzupełnienie macierzy A o wszystkie wektory bi:

26 Dekompozycja LU W fazie eliminacji Gaussa-Jordana przekształceniom podlegają wszystkie wektory bi, Faza wstecznego podstawiania jest realizowana dla każdego bi osobno Dekompozycja LU oferuje bardziej ekonomiczny i elastyczny schemat (nie trzeba znać z góry bi)

27 Dekompozycja LU Dekompozycja LU kwadratowej macierzy A: gdzie L – macierz trójkątna dolna, U – macierz trójkątna górna:

28 Dekompozycja LU Dekompozycja nie jest unikalna (nieskończenie wiele możliwych) Niektóre ograniczenia dają metodę z nazwą:

29 Dekompozycja LU Idea zastosowania:Najpierw rozwiązanie Ly=b (podstawianie wprzód) Następnie rozwiązanie Ux=y (podstawianie wstecz) Gdy już są znane L, U, to dla każdego wektora bi potrzebne są dwa tanie procesy podstawiania

30 Dekompozycja DoolittleDekompozycja ma ścisły związek z eliminacją Gaussa Po dekompozycji: Przed dekompozycją:

31 Dekompozycja DoolittleZastosowanie eliminacji Gaussa: A(2,:)=A(2,:)-L21 A(1,:) A(3,:)=A(3,:)-L31 A(1,:) A(3,:)=A(3,:)-L32 A(2,:)

32 Dekompozycja DoolittleMacierz U uzyskiwana w trakcie eliminacji Gaussa-Jordana jest identyczna, jak w dekompozycji Doolittle Macierz L dekompozycji Doolittle składa się ze współczynników używanych do eliminacji kolejnych równań w metodzie Gaussa-Jordana Dekompozycję Doolittle można więc zrealizować poprzez eliminację Gaussa-Jordana, w której mnożniki używane do eliminacji tworzą macierz L

33 Dekompozycja DoolittleMnożniki eliminacyjne można zapisywać w miejsce zer w przekształcanej macierzy A

34 Dekompozycja Doolittle

35 LU: rozwiązanie układu I

36 LU: rozwiązanie układu II

37 Dekompozycja Cholesky’egoDekompozycja o postaci: Macierz A musi być: Symetryczna Dodatnio określona Takie macierze pojawiają się w wielu zastosowaniach, jak choćby w metodzie najmniejszych kwadratów (później w tym semestrze)

38 Cholesky

39 Cholesky Macierz A jest symetryczna, więc można się zajmować tylko częścią trójkątną dolną Sześć równań

40 Cholesky pierwsza kolumna: druga kolumna: trzecia kolumna:

41 Cholesky Ogólniej:

42 Cholesky Dla pierwszej kolumny: Dowolna kolumna – niewiadomą jest Lij:Dla elementu diagonalnego (i=j):

43 Cholesky

44 Specjalne rodzaje macierzyMacierze rzadkie – większość elementów ma wartość zero Macierze wstęgowe – macierze rzadkie, w których elementy niezerowe są położone tylko na głównej przekątnej i przekątnych przyległych

45 Macierze rzadkie Macierz trójprzekątniowa – elementy tylko na głównej przekątnej i dwóch przekątnych przyległych

46 Macierz trójprzekątniowaBardzo mało kosztowne pod względem pamięci (np. dla m=n=100, 298 vs liczb)

47 Dekompozycja LU 3-diag A(k,:)=A(k,:)-c(k-1)/d(k-1)*A(k,:), k=2,3...ne(k) nie ulega zmianie, mnożnik eliminacyjny jest zapisywany w pozycji c(k-1):

48 Dekompozycja LU 3-diag

49 Rozwiązanie po dekompozycjiFaza Ly=b:

50 Rozwiązanie po dekompozycjiFaza Ux=y:

51 Metody iteracyjne

52 Algorytm Gaussa-SeidelaAlgebraiczny zapis i-tego równania Wyłączając komponent xi:

53 Algorytm Gaussa-SeidelaRozwiązując: Pomysł na schemat iteracyjny:

54 Algorytm Gaussa-SeidelaProcedura startuje z pewnego wstępnie wybranego (np. losowo) oszacowania rozwiązania xo Zastosowanie wzoru: poprawia oszacowanie i tworzy jeden cykl algorytmu Procedura jest powtarzana dotąd, aż zmiana xp w stosunku do xp-1 w kolejnej iteracji p jest nieznaczna

55 Gauss-Seidel - relaksacjaWartość xp w kolejnej iteracji jest sumą ważoną: ω jest współczynnikiem relaksacji: Gdy ω=1 – efektywnie brak relaksacji Gdy ω < 1 – interpolacja Gdy ω >= 1 - ekstrapolacja

56 Gauss-Seidel

57 Dobieranie współczynnika ωZazwyczaj stosuje się przepis: gdzie: jest przeregulowaniem w k-tej iteracji bez relaksacji, k powinno być rzędu 10, a p >= 1

58 Gauss-Seidel z relaksacjąPrzeprowadzić k (np. k=10) iteracji bez relaksacji, czyli dla ω = 1, zarejestrować Δx(k) Przeprowadzić dodatkowe p (p >= 1) iteracji bez relaksacji, zarejestrować Δx(k+p) Z podanego wzoru wyznaczyć ωopt i resztę obliczeń przeprowadzać z relaksacją ωopt

59 Piwot - motywacja Kolejność używania równań może mieć decydujący wpływ na rozwiązanie

60 Piwot - motywacja Przykład. (rozwiązanie  x1=x2=x3=1„Dobre” uporządkowanie: Algorytm Gaussa-Jordana zawiedzie (dzielenie przez zero w pierwszej operacji) „Złe” uporządkowanie:

61 Dwa ostatnie równania są sprzecznePiwot - motywacja A gdyby to nie było zero, ale bardzo mała wartość... 1/ε, 2/ε jest bardzo duże Dwa ostatnie równania są sprzeczne Po eliminacji pierwszego równania:

62 Piwot - motywacja Jeśli nawet zła kolejność rozwiązywania nie prowadzi do takich krytycznych sytuacji, to akumulacja błędów zaokrąglania czyni często rozwiązanie całkowicie bezwartościowym

63 Diagonalna dominacja Stosując piwot, dążymy do uzyskania diagonalnej dominacji Macierz wykazuje d.d., jeśli w każdym wierszu wyraz na przekątnej jest większy od sumy pozostałych wyrazów (w sensie wartości bezwzględnych):

64 Diagonalna dominacja Przykład: Macierz nie wykazuje d.d.:... ale po przestawieniu wierszy wykazuje:

65 Skalowany piwot Jak poprzestawiać wiersze, żeby uzyskać d.d. (wtedy najmniej podatne na problemy rozwiązanie) Przydatne też, aby wyraz na przekątnej miał nie tylko wartość dominującą w wierszu, ale i możliwie jak największą z wszystkich dostępnych w macierzy

66 Skalowany piwot Dodatkowy wektor kolumnowy współczynników skalujących:Wyznaczony przez:

67 Skalowany piwot Rozstrzygnięcia posługują się relatywną wielkością wyrazów: Eliminacja przy użyciu wiersza k: Zamiast akceptować automatycznie Akk jako dzielnik, poszukiwany jest poniżej lepszy wiersz piwotowy (dlaczego tylko poniżej?)

68 Skalowany piwot Wybierany jest wiersz p, w którym Apk ma największą względną wielkość: Po znalezieniu takiego wiersza, zamienia się go miejscami z wierszem k (stosowna zamiana musi też być dokonana w wektorze współczynników skalujących s)

69 Skalowany piwot

70 Skalowany piwot Pomocnicza funkcja swapRows:

71 Dekompozycja DoolittleW piwocie zamieniamy miejscami wiersze, dbając to, żeby zamiana dotyczyła także wektora b (RHS) Przy dekompozycji LU szykujemy się na „przyszłe” wektory bj, ale ich wyrazy będą uporządkowane wg pierwotnej kolejności Wyrazy wektorów bj przed rozwiązaniem korzystającym z LU należy poddać przestawieniu takiemu samemu, jak w czasie eliminacji Gaussa

72 Dekompozycja DoolittleW algorytmie Gaussa z piwotem potrzebny jest dodatkowy wektor permutacyjny p, wstępnie równy (1:n)’ Zamiana wierszy miejscami dotyczy także wektora p Po zakończonym algorytmie G.-J. wektor p wskazuje kolejność wierszy Uporządkowanie „przyszłych” wektorów b  w zapisie Matlab: b(p)