Rozwiązywanie układów równań liniowych

Rozwiązywanie układów równań liniowych

Szukamy rozwiązania układu równań liniowych w postaci:

Powyższy układ rownań można zapisać w postaci macierzowej:

gdzie:

- szukany

- macierz

współczynników układu

równań

- wektor

Rozwiązywanie układów równań liniowych

Układ równań może mieć nieskończenie wiele rozwiązań, dokładnie jedno, lub rozwiązanie może nie istnieć wcale.

Pojęcia podstawowe.

Def. Jeżeli macierz A jest macierzą kwadratową to macierz B jest macierzą do niej odwrotną jeśli zachodzi pomiędzy nimi związek:

Def. Macierz jednostkowa jest macierzą kwadratową, której elementy leżące na diagonali są równe 1 a pozostałe 0.

Jeżeli nie istenieje taka macierz B to macierz A jest nieodwracalna.

Macierz kwadratowa o wyznaczniku różnym od 0 jest odwracalna i nazywamy ją nieosobliwą. Macierz osobliwa ma wyznacznik równy 0.

Rozwiązywanie układów równań liniowych

Def. Transpozycję macierzy A dokonuje się poprzez zamianę wierszy w kolumny a kolumn w wiersze. Macierz transponowaną oznaczamy A^T

Własności macierzy transponowanej:

Transpozycja macierzy nie wpływa na wyznacznik i na ślad macierzy:

Jeżeli a i b są wektorami to zachodzi relacja:

Rozwiązywanie układów równań liniowych

Def. Macierz kwadratową A nazywamy idempotentną jeśli zachodzi związek:

Def. Macierzą dodatnio określoną nazywamy macierz o własności:

1) Jeśli A jest macierzą zespoloną hermitowską (a_ij=a_ji^*) to dla x∈Rⁿ :

2) Jeśli macierz A jest macierzą rzeczywistą to dla x∈Rⁿ :

Macierz dodatnio określona jest zawsze odwracalna. Macierz odwrotna jest również dodatnio określona.

Mnożenie macierzy nie jest w ogólności przemienne:

Rozwiązywanie układów równań liniowych

Dla macierzy A stopnia n ślad macierzy definiujemy następująco:

Normy wektorów i macierzy Dla n-wymiarowego wektora

najczęściej stosowane są poniższe normy:

1) 2)

3)

Dla dowolnego wektora x w przestrzeni Rⁿ prawdziwe są poniższe relacje pomiędzy normami:

Rozwiązywanie układów równań liniowych

Macierz o m wierszach i n kolumnach można traktować jako operator liniowy przekształcający przestrzeń R^m w Rⁿ. Normę takiej macierzy można określić przy użyciu wektorów:

gdzie: p i q oznaczają normy wektorów w przestrzeniach Rⁿ i w R^m

jest normą macierzy indukowaną przez normy wektorów Dla p=q oznaczając:

możemy określić następujące normy macierzy:

- maksymalna suma modułów w kolumnie

- największa wartość własna macierzy AA^T

- maksymalna suma modułów

Rozwiązywanie układów równań liniowych

- maksymalny moduł elementu

W przestrzeniach z normą |⋅|² często używa się euklidesowej normy macierzy:

Nie jest ona indukowana żadną normą ale spełnia ona z normą ||⋅||² warunek zgodności:

Normy macierzy mają istotne znaczenie w analizie błędów rozwiązania (np. układów równań liniowych).

Rozwiązywanie układów równań liniowych

Układ równań z macierzą trójkątną.

Zakładamy, że elementy leżące na diagonali są niezerowe. Rozwiązanie układu można znaleźć posługując się wzorem rekurencyjnym, zaczynając od elementu x_n:

Najprostszy (najłatwiejszy do rozwiązania) układ równań ma postać:

W celu wyznacznenia wszystkich składowych wektora rozwiązania x należy wykonać:

- operacji mnożenia i dzielenia

Rozwiązywanie układów równań liniowych

Podczas wykonywania obliczeń pojawiają się błędy zaokrągleń wpływające na wynik.

Błędy te są równe zastępczemu zaburzeniu , które spełnia nierówność:

Z oszacowania tego wynika, że:

Znając oszacowanie od góry elementów macierzy A, tj.

można oszacować prawe strony nierówności:

Rozwiązywanie układów równań liniowych

W przypadku gdy macierz jest diagonalnie dominująca (wiersz) tzn. spełniony jest warunek:

wówczas oszacowanie jest następujące:

Natomiast gdy macierz współczynników układu jest

diagonalnie dominująca kolumnowo, otrzymujemy oszacowanie:

Rozwiązywanie układów równań liniowych

Metoda eliminacji Gaussa rozwiązania układu równań liniowych.

Metoda jest dwuetapowa:

1) Eliminacja zmiennych.

Wyjściowe równanie ma postać:

Odejmujemy od i-tego wiersza (i=2,3,...,n) wiersz pierwszy pomnożony przez współczynnik

Z równań i=2,3,..,n wyeliminowana została zmienna x .

Rozwiązywanie układów równań liniowych

Powtarzamy operację, ale odejmujemy od i-tego wiersza (i=3,4,...,n) wiersz drugi pomnożony przez współczynnik

Postępując dalej w ten sposób eliminujemy z każdego następnego równania jedną zmienną. Eliminację kończymy po (n-1) krokach, gdy uzyskamy trójkątny układ równań w postaci:

2) Etap drugi nazywany jest postępowaniem odwrotnym. Rozwiązanie (kolejne składowe wektora x) znajdujemy stosując wzór rekurencyjny dla macierzy trójkątnej.

Wyznaczenie rozwiązania metodą Gaussa wymaga wykonania:

- operacji mnożenia i dzielenia

- operacji dodawania i odejmowania

Rozwiązywanie układów równań liniowych

Uwaga: Podczas obliczeń może wystąpić problem dzielenia przez zero (gdy element a_ii=0). W tkim przypadku należy zmienić kolejność wierszy, najlepiej

tak aby dzielić przez element o największym module. Operacja zmiany kolejności wierszy nosi nazwę wyboru elementu podstawowego. Należy jednak pamiętać żeby zmienić także kolejność elementów w wektorze wyrazów wolnych oraz

w wektorze rozwiązań x.

Metoda eliminacji Jordana (eliminacji zupełnej) W układzie równań:

równanie pierwsze dzielimy obustronnie przez współczynnik:

i odejmujemy od i-tego wiersza (i=2,3,...,n) wiersz pierwszy przemnożony przez

Rozwiązywanie układów równań liniowych

Otrzymujemy wówczas układ równań:

Podobnie postępujemy z równaniem drugim. Dzielimy je obustronnie przez:

Następnie od i-tego wiersza (i=1,3,4,...,n) odejmujemy wiersz drugi pomnożony przez współczynnik:

Otrzymujemy zmodyfikowany układ równań:

Rozwiązywanie układów równań liniowych

Po przeprowadzeniu (n-1) eliminacji zmiennych układ rówań ma poniższą postać:

czyli otrzymujemy gotowe rozwiązanie.

Metoda eleminacji Jordana wymaga wykonania:

Nakład obliczeń jest 1.5 razy większy niż w metodzie eliminacji Gaussa.

Jeśli pojawia się problem dzielenia przez 0, należy dokonać wyboru elementu podstawowego.

Rozwiązywanie układów równań liniowych

Rozkład LU metodą Gaussa(-Crouta)

Metodę Gaussa można użyć do znalezienia takch macierzy L i U, które z macierzą A związane są relacją:

Procedura wyznaczania elementów tych macierzy nosi nazwę rozkładu LU.

Sposób postępowania (wykorzystujemy metodę eliminacji Gaussa):

1) mnożenie wierszy (i=2,3,...,n) w układzie równań przez czynnik

i odjęcie od i-tego wiersza, wiersza o indeksie 1 zastępuje się mnożeniem przez macierz:

Co można zapisać macierzowo:

Rozwiązywanie układów równań liniowych

Eliminacja zmiennej z równań (i=3,4,...,n) wygląda podobnie. Mnożymy wiersze zmodyfikowanego układu równań o indeksach i=3,4,...,n przez czynnik

i odejmujemy od nich wiersz drugi. Operację tą można przeprowadzić mnożąc układ równań obustronnie przez macierz L⁽²⁾:

Zapis macierzowy operacji:

Po wykonaniu (n-1) razy takich operacji otrzymamy:

Rozwiązywanie układów równań liniowych

Macierze L⁽ⁱ⁾ są nieosobliwe (można znaleźć dla każdej macierz odwrotną).

Przemnażając obie strony powyższych równań przez (L^(n-1))^-1, (L^(n-2))^-1...., otrzymamy:

Wprowadzamy oznaczenia:

Należy jeszcze znaleźć macierze odwrotne (L⁽ⁱ⁾)^-1:

Rozwiązywanie układów równań liniowych

Macierz L:

Uwaga: Macierz L jest macierzą dolną z jedynkami na diagonali, macierz U jest macierzą trójkątną górną.

Dysponując macierzami L i U można rozwiązać układ równań:

Rozwiązujemy najpierw szukając wektora y:

Rozwiązywanie układów równań liniowych

Rozkład LU metodą Doolittle'a.

Równanie A=LU traktujemy jako n² równań z n² niewiadomymi l_ij (i>j) i niewiadomymi u_ij (i≤j ). Elementy l_ij oraz u_ij oblicza się ze wzorów:

Liczba wykonywanych działań jest identyczna jak w metodzie Gaussa.

Elementy u_ij oraz l_ij można zapamiętać w komórkach macierzy A.

Rozwiązywanie układów równań liniowych

Zalety:

1) Duża wydajność dla dużej liczby równań.

Rozkład LU opłaca się stosować w przypadku rozwiązywania wielu układów równań z tą samą macierzą współczynników układu A. Każdy układ równań różni się wtedy tylko wektorem wyrazów wolnych. Rozkład LU wykonuje się w takim przypadku tylko raz (ilość operacji ~n³). Rozwiązanie pojedynczego układu równań można znaleźć przy zastosowaniu algorytmu postępowania odwrotnego (ilość operacji ~n²).

2) Oszczędność zajmowanej pamięci.

Elementy macierzy L i U mogą zostać zapisane w macierzy A. Pomija się elementy diagonalne, których wartośc jest równa 1.

3) Jeśli macierz A jest symetryczna i dodatniookreslona to nie trzeba dokonywać wyboru elementów podstawowych. Jeśli macierz jest symetryczna i silnie

diagonalnie dominująca to macierz jest dodatniookreślona.

Rozwiązywanie układów równań liniowych

Układy równań z macierzą symetryczną.

Rozkład LDL^T

Oznaczmy rozkład LU jako:

Szukamy rozkładu macierzy A w postaci:

gdzie: L – macierz trójkątna dolna z jedynkami na diagonali

D – macierz diagonalna z elementami diagonalnymi macierzy U – macierz trójkątna górna z jedynkami na diagonali

Jeśli macierz A jest symetryczna:

Dla macierzy symetrycznych istnieje jednoznaczny rozkład:

Rozwiązywanie układów równań liniowych

Elementy macierzy L i D wyznaczamy ze wzorów:

a dla i=2,3,...,n oblicza się na przemian:

W celu wyznaczenia rozkładu LDL^T należy wykonać:

-operacji mnożenia i dzielenia

-operacji dodawania i odejmowania

Rozwiązywanie układów równań liniowych

Układ równań Ax=b rozwiązuje się w trzech etapach:

Zaletą rozkładu LDL^T jest dwukrotnie mniejszy nakład obliczeń potrzebny do znalezienia rozwiązania układu równań niż w metodzie Gaussa-Crouta.

Ponadto zamiast posługiwać się pełną macierzą o n² elementach wystarczy zapamiętać dolną macierz trójkątną z elementami diagonalnymi. Liczba elementów potrzebnych do zapamiętania wynosi zatem:

Rozwiązywanie układów równań liniowych

Rozkład Banachiewicza-Cholesky'ego

Jeśli macierz A jest macierzą symetryczną dodatnio okresloną wówaczas można dokonać następującego rozkładu:

Macierz L jest macierzą trójkątną dolną z elementami na diagonali mogącymi się różnić od 1.

Macierz

z rozkładu LU również spełnia warunek:

więc rozkład ten nie jest jednoznaczny. Jeśli jednak liczby na diagonali macierzy L są ddatnie wówczas rozkład jest jednoznaczny, a elementy macierzy wyznaczamy ze wzorów:

Rozwiązywanie układów równań liniowych

Aby przeprowadzić rozkład LL^T należy wykonać:

n obliczeń pierwiastka kwadratowego

- operacji mnożenia i dzielenia

- operacji dodawania i odejmowania Rozkład LL^T pozwala na wykonanie dwukrotnie mniejszej liczby operacji oraz zajęcie prawie dwukrotnie mniejjszej liczby komórek pamięci niż metoda Gaussa-Crouta.

Rozwiązywanie układów równań liniowych

Przykład:

obliczenia elementów:

Rozwiązywanie układów równań liniowych

Układy równań z macierzą trójdiagonalną Szukamy rozwiązania układu równań:

Zdarza się że macierz układu równań jest macierzą postaci:

Można wykonać rozkład LU macierzy T, macierze te mają postać:

Rozwiązywanie układów równań liniowych

Elementy l_i oraz u_i oblicza się rekurencyjnie:

Rozwiązanie układu Tx=d czyli Lux=d rozwiązuje się dwuetapowo:

wykorzystując wzory:

Rozwiązywanie układów równań liniowych

Dokonując rozkładu T=LU należy wykonać:

M=2n-2 - operacji mnożenia i dzielenia oraz

D=n-1 - operacji dodawania i odejmowania

Uwzględniając dodatkowo rozwiązanie układu równań należy w sumie wykonać:

M=3n-2 i D=2n-2 działań. Ponadto elementy macierzy L i U można zapamiętać

w macierzy T. Liczba zajętych komórek pamięci wynosi: P=3n-2. Zamiast posługiwać się macierzami można używać 3 wektorów jednowymiarowych.

Jeśli macierz jest dominująca kolumnowo to rozkład T=LU jest równoważny

rozkładowi z częściowym wyborem elementu podstawowego (niezawodność metody).

Elementy macierzy L spełniają warunek:

co poprawia oszacowanie błędów zaokrągleń.

Rozwiązywanie układów równań liniowych

Rozkład LU – obliczanie wyznacznika

Aby obliczyć wyznacznik macierzy A możemy posłużyć się rozkładem

Macierz E jest macierzą zaburzeń. Jej obecność w powyższych wzorach jest wynikiem błędów zaokrągleń pojawiających się podczas obliczania rozkładu.

Wyznacznik macierzy U można łatwo obliczyć, gdzyż jest on iloczynem n elementów stojących na diagonali tej macierzy. Dysponując rozkładem LU wystarczy wykonanie n-1 operacji mnożenia w celu wyznaczenia wyznacznika macierzy A.

Małe zaburzenia wartości elementów macierzy mogą powodować duże zmiany wartości wyznacznika macierzy. Zmiany te zależą od wielkości:

Rozwiązywanie układów równań liniowych

Rozkład LU – odwracanie macierzy.

Aby znaleźć przy pomocy macierzy L i U macierz odwrotną A^-1 należy rozwiązać n układów równań:

gdzie: -jest wersorem w n-wymiarowej przestrzeni

Rozwiązania układów równań x⁽ⁱ⁾ stanowią kolumny macierzy odwrotnej A^-1 (po uwzględnieniu ewentualnych przestawień wierszy wynikających z

wyboru elementu podstawowego).

Przykład. Znaleźć macierz A^-1 jeśli macierz A jest zdefiniowana:

Rozwiązywanie układów równań liniowych

Wektor przestawień (wektor permutacji)

Macierz odwrotna:

Rozwiązywanie układów równań liniowych

Iteracyjne poprawianie rozwiązania układu równań Błąd rozwiązania można sprawdzić obliczając wektor reszt:

Zazwyczaj współrzędne wektora r są różne od zera. Oznacza to że nie uzyskaliśmy dokładnego rozwiązania, ale przybliżone. Rozwiązanie to chcemy poprawić:

gdzie: jest poprawką, którą można łatwo wyznaczyć

Należy jednak pamiętać że wyznaczając poprawkę znowu popełniamy jakiś błąd i rozwiązanie które uzyskamy będzie miało postać:

Rozwiązywanie układów równań liniowych

Lemat. Jeżeli wektor reszt r=b-Ax jest obliczony dokładnie, poprawka x została wyznaczona metodą Gaussa-Crouta oraz zachodzi warunek:

to wtedy:

Sposób iteracyjnego poprawiania rozwiązania:

1. Rzowiązujemy układ Ax⁽¹⁾=b metodą Gaussa

2. Obliczamy wektor reszt r⁽¹⁾ i sprawdzamy rozwiązanie 3. Sprawdzamy czy poniższy warunek jest prawdziwy

jeśli nie to przerywamy obliczenia. Jeśli jest spełniony to kontynuujemy.

4. Obliczamy poprawkę i wyznaczamy

5. Wyznaczamy wektor reszt r⁽²⁾ i sprawdzamy rozwiązanie. W razie konieczności

Rozwiązywanie układów równań liniowych

Metody iteracyjne rozwiązywania układów równań.

Idea metody. Rozważamy ciąg wektorów x⁽⁰⁾, x⁽¹⁾,x⁽⁰⁾,...generowany przez wyrażenie rekurencyjne:

gdzie: x⁽⁰⁾ jest ustalonym wektorem, M jest macierzą kwadratową, w jest wektorem.

Ciąg określony powyższą relacją rekurencyjną przy dowolnym wektorze x⁽⁰⁾ jest zbieżny do jedynego punktu granicznego wtedy i tylko wtedy gdy:

Promień spektralny jest to największa co do wartości bezwględnej wartość własna macierzy.

Poszukując iteracyjnie kolejnych elementów ciągu dążymy do osiągnięcia:

Aby znaleźć iteracyjnie rozwiązanie układu równań należałoby użyć dowolnej macierzy M okreslonej w lemacie oraz wektor w zdefiniowany:

Rozwiązywanie układów równań liniowych

Takie postępowanie wymagałoby znalezienia macierzy odwrotnej układu co byłoby nieekonomiczne. Aby tego uniknąć zakładamy:

gdzie: N jest pewną macierzą kwadratową.

Określamy macierz M w postaci:

i otrzymujemy ogólny przepis iteracyjny:

Zastosowanie powyższej realacji pozwala wyznaczyć wektor rozwiązań jeśli spełniony jest warunek:

Istnieją sposoby określenia macierzy N. istotnym problemem jest efektywność poszczególnych metod (liczba wykonywanych działań ~ liczba iteracji)

oraz błędy zaokrągleń wpływające na wynik końcowy. Efektywność można badać określając wektor błędu:

Rozwiązywanie układów równań liniowych

Efektywność metod iteracyjnych związana jest ściśle z promieniem spektralnym macierzy M.

Metody iteracyjne są bardzo czułe na błędy zaokrągleń. W każdej iteracji zamiast:

otzrymujemy:

Nawet przy małej wartości błędów zaokrągleń, duża liczba iteracji może spowodować kumulację błędów:

W skrajnym przypadku błąd zaokrągleń:

może doprowadzić do uzyskania

i ciąg wektorów nie jest zbieżny do rozwiązania.

Rozwiązywanie układów równań liniowych

Metoda Jacobiego.

W metodzie jacobiego macierz układu A rozkładamy na 3 składniki

gdzie: L – macierz poddiagonalna, D – macierz diagonalna, U – macierz naddiagonalna

Jeśli przyjmiemy:

wówczas macierz M ma postać:

Rozwiązywanie układów równań liniowych

Wzór iteracyjny Jacobiego przyjmuje postać:

Elementy D^-1 są odrotnościami elementów diagonalnych macierzy D.

Aby móc zastosować wzór Jacobiego należy najpierw sprawdzić czy w macierzy A na diagonali nie ma elementów zerowych. Jeśli są, wówczas należy przestawić wiersze w A tak by miały one niezerowe wartości. Procedura zmiany kolejności wierszy nie gwarantuje spełnienia waruku zbieżności:

Metoda jest zbieżna jeśli macierz jest silnie diagonalnie dominująca lub silnie diagonalnie dominująca kolumnowo.

Rozwiązywanie układów równań liniowych

Metoda Gaussa-Seidla

Zakładamy że macierz układu jest sumą 3 macierzy:

ale przyjmujemy inną postać macierzy N:

Macierz M w procedurze iteracyjnej wyraża się następująco:

Wzór iteracyjny w metodzie Gaussa-Seidla:

Rozwiązywanie układów równań liniowych

Nie ma problemu z obliczeniem prawej strony równania. Zaczynając obliczenia

od elementu pierwszego tj. x₁⁽ⁱ⁺¹⁾ , można zauważyć, że obliczana wartość nie zależy

od x⁽ⁱ⁺¹⁾ (zapewnia to macierz L). Przy wyznaczaniu wartości x₂⁽ⁱ⁺¹⁾ zauważamy, że

potrzebujemy x₁⁽ⁱ⁺¹⁾ itd. Obliczając i-ty element wektora x⁽ⁱ⁺¹⁾ korzystamy

z elementów już obliczonych czyli, x_k⁽ⁱ⁺¹⁾, i=1,2,...,i-1. Zabieg ten może w pewnych przypadkach znacznie podnieść efektywność metody.

Metodę GS można zastosować w przypadku niezerowych elementów diagonalnych macierzy A. Metoda jest zbieżna jeśli macierz jest symetryczna i dodatnio określona oraz gdy jest silnie diagonalnie dominująca lub silnie diagonalnie dominująca

kolumnowo. Na ogół metoda Gaussa-Seidla jest szybciej zbieżna od metody Jacobiego, choć nie jest to regułą.

Nakład obliczeń i testy stopu.

Zadaniem metod iteracyjnych jest zmniejszenie początkowego błędu do małej wartości 

Nie można wieć z góry przewidzieć jaka będzie potrzebna liczba iteracji aby spełnić ten warunek. Liczbę ta można jedynie oszacować:

Rozwiązywanie układów równań liniowych

Oszacowanie to nie uwzględnia błędów zaokrągleń, które mogą zmniejszyć efektywność metody.

W celu uniknięcia niedoszacowania rozwiązania należy zastosować test stopu:

a)

b)

Testy te mogą być jednak zawodne. Jeśli norma macierzy A jest mała to wartość reszty

Może być mała mimo dużego odchylenia x⁽ⁱ⁺¹⁾ od rozwiązania dokładnego.

Również mała wartość

nie świadczy o zbliżaniu się do rozwiązania dokładnego, co wynika z poniższego

Rozwiązywanie układów równań liniowych

Jeśli macierz (I-M) ma małą normę to mimo dużego błędu e⁽ⁱ⁾ wektory

x⁽ⁱ⁺¹⁾ oraz x⁽ⁱ⁾ mogą się niewiele różnić. W przypadku wolnej zbieżności

należy wówczas zastosować test:

lub jeśli potrafimy oszacować normę macierzy A^-1 test: