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:

W dokumencie Rozwiązywanie układów równań liniowych (Stron 21-44)