Metody numeryczne fizyki/Algebraiczne sposoby rozwiązywania układów równań liniowych

Algebraiczne sposoby rozwiązywania układów równań liniowych

---

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, na tych samych warunkach, z możliwością obowiązywania dodatkowych ograniczeń. 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ść.

---

---

Będziemy się tutaj zajmowali rozwiązaniem algebraicznych układów równań, które możemy przedstawić w postaci macierzowej $\mathbf {A} \mathbf {x} =\mathbf {b} \;$, gdzie $\mathbf {A} \;$ jest macierzą o m wierszach i n kolumnach, a $\mathbf {x} \;$ jest wektorem o n-niewiadomych, i $\mathbf {b} \;$ jest wektorem o n wyrazach wolnych. Jeśli z tego równania macierzowego wyznaczymy wektor niewiadomych, to możemy wyznaczyć ten wektor znając macierz $\mathbf {A} \;$ i macierz wyrazów wolnych. Jest to ścisły sposób wyznaczania wektora niewiadomych. Istnieją też numeryczne metody takiego wyznaczania wektora niewiadomych, które poniżej przedstawimy.

Wprowadzenie do pojęcia normy

Weźmy sobie przestrzeń $\mathbf {R} ^{n}\;$, którego elementami są wektory pionowe $\mathbf {x} =[x_{1},x_{2},..,x_{n}]^{T}\;$, wtedy możemy wprowadzić poszczególne normy inaczej zdefiniowane stosowanych w obliczeniach numerycznych:

$||\mathbf {x} ||_{1}=|x_{1}|+|x_{2}|+..+|x_{n}|\;$

(5.1)

$||\mathbf {x} ||_{2}=\left(x_{1}^{2}+x_{2}^{2}+...+x_{n}^{2}\right)^{1/2}\;$

(5.2)

$||\mathbf {x} ||_{\infty }=\max\{|x_{1}|,|x_{2}|,..,|x_{n}|\}\;$

(5.3)

Normy (5.1), (5.2) i (5.3) spełniają warunki poniżej dla dowolnego wektora należącego do n-wymiarowej przestrzeni rzeczywistej, które to piszemy:

$||\mathbf {x} ||_{\infty }\leq ||\mathbf {x} ||_{2}\leq ||\mathbf {x} ||_{1}\leq {\sqrt {n}}||\mathbf {x} ||_{2}\leq n||\mathbf {x} ||_{\infty }\;$

(5.4)

Określmy teraz macierz o m wierszach i n kolumnach, która może być traktowana jako operator przekształcający wektor należącej do n-wymiarowej przestrzeni rzeczywistej na m-wymiarową przestrzeń rzeczywistą, wtedy możemy określić normę tej naszej omawianej macierzy $\mathbf {A} \;$:

$||\mathbf {A} ||_{pq}=\max _{\overset {x\in \mathbf {R} ^{n}}{\underset {x\neq 0}{}}}{{||\mathbf {A} \mathbf {x} ||_{q}} \over {||\mathbf {x} ||_{p}}}\;$

(5.5)

Symbole "p" i "q" oznaczają normy przedstawione w punktach (5.1), (5.2) i (5.3), dla którego definicja normy macierzy jest przedstawiona w punkcie (5.5). Zdefiniujmy teraz trzy kolejne normy macierzy, podobne do definicji norm wektora $\mathbf {x} \;$ (5.1), (5.2) i (5.3), tylko że ich definicje wyglądają:

$||\mathbf {A} ||_{1}=\max _{j=1,2,..,n}\sum _{i=1}^{m}|a_{ij}|\;$ (5.6)

$||\mathbf {A} ||_{\infty }=\max _{i=1,2,..,m}\sum _{j=1}^{n}|a_{ij}|\;$ (5.7)

$||\mathbf {A} ||_{1\infty }=\max _{ij}|a_{ij}|\;$

(5.8)

A także zdefiniujmy $||\mathbf {A} ||_{2}\;$ jako największa wartość własną macierzy $(\mathbf {A} ^{T}\mathbf {A} )^{1/2}\;$. Wszystkie powyższe normy operatora $\mathbf {A} \;$ możemy obliczyć łatwo numerycznie, tylko $||\mathbf {A} ||_{2}\;$ jest trudne do obliczenia. Przestrzeń z normą ||⋅||₂ używa sie go w eulkidesowej normy macierzy, które często są zwane normami Schura lub normami Frobeniusza, którego definicja jest:

$||\mathbf {A} ||_{E}={\sqrt {\sum _{i=1}^{m}\sum _{j=1}^{n}a_{ij}^{2}}}\;$

(5.9)

wtedy możemy napisać warunek zgodności dla n-wymiarowej przestrzeni rzeczywistej, do których należą elementy $\mathbf {x} \;$, zatem:

$||\mathbf {A} \mathbf {x} ||_{2}\leq ||\mathbf {A} ||_{E}||\mathbf {x} ||_{2}\;$

(5.10)

Dla dowolnej macierzy $\mathbf {A} \;$ możemy powiedzieć, że zachodzi warunek, który jest zawsze prawdziwy:

$||\mathbf {A} ||_{2}\leq ||A||_{E}\leq {\sqrt {n}}||A||_{2}\;$

(5.11)

Dla dowolnych dwóch macierzy i dla dowolnej definicji normy możemy napisać twierdzenie dla dowolnych indeksów p,q,r, która jest również słuszna dla dowolnej normy:

$||\mathbf {A} \mathbf {B} ||_{pq}=\max _{\overset {x\in \mathbf {R} ^{n}}{\underset {x\neq 0}{}}}{{||\mathbf {A} \mathbf {B} \mathbf {x} ||_{q}} \over {||\mathbf {x} ||_{p}}}=\max _{\overset {x\in \mathbf {R} ^{n}}{\underset {x\neq 0}{}}}{{||\mathbf {A} \mathbf {B} \mathbf {x} ||_{q}} \over {||\mathbf {B} \mathbf {x} ||_{r}}}{{||\mathbf {B} \mathbf {x} ||_{r}} \over {||\mathbf {x} ||_{p}}}\leq \max _{\overset {x\in \mathbf {R} ^{n}}{\underset {x\neq 0}{}}}{{||\mathbf {A} \mathbf {B} \mathbf {x} ||_{q}} \over {||\mathbf {B} \mathbf {x} ||_{r}}}\max _{\overset {x\in \mathbf {R} ^{n}}{\underset {x\neq 0}{}}}{{||\mathbf {B} \mathbf {r} ||_{q}} \over {||\mathbf {x} ||_{p}}}=\;$
$=||\mathbf {A} ||_{rq}||\mathbf {B} ||_{pr}\Rightarrow ||\mathbf {A} \mathbf {B} ||_{pq}\leq ||\mathbf {A} ||_{rq}||\mathbf {B} ||_{pr}\;$

(5.12)

Również możemy powiedzieć $||A||_{pr}\leq ||A^{'}||_{pr}\;$ z definicji macierzy podobnych, że normy macierzy podobnych mają równe wartości. Można też udowodnić, że jeśli λ_i jest jedną z wartości własnych macierzy $\mathbf {A} \;$, która może być również wartością zespoloną przy zdefiniowanej macierzy kwadratowej z pewną normą wektora nałożoną na macierz $\mathbf {A} \;$, to możemy napisać:

$||\mathbf {A} \mathbf {x} ||=|\lambda _{i}|||\mathbf {x} ||_{p}\Rightarrow ||\mathbf {A} ||_{pr}||x||\geq |\lambda _{i}||x||_{p}\Rightarrow |\lambda _{i}|\leq ||\mathbf {A} ||_{pr}\;$

(5.13)

Wykorzystując definicję normy $||\cdot ||_{2}\;$, możemy napisać wniosek przy pomocy normy zdefiniowanej na macierzy $\mathbf {A} \;$, wykorzystując twierdzenia (5.12):

$||\mathbf {A} ||_{2}^{2}\leq ||\mathbf {A} ^{T}\mathbf {A} ||_{\infty }\leq ||\mathbf {A} ^{T}||_{\infty }||\mathbf {A} ||_{\infty }=||\mathbf {A} ||_{1}||\mathbf {A} ||_{\infty }\;$

(5.14)

Zdefiniujemy teraz dwa lematy, który są dla nas ważne wykorzystując niektóre dowodu napisane powyżej:

Lemat pierwszy
Jeśli norma macierzy $\mathbf {M} \;$ spełnia warunek $||\mathbf {M} _{n\times n}||<1\;$, to macierz $\mathbf {I} +\mathbf {M} \;$ jest macierzą nieosobliwą spełniającej warunek dla p=1,2,∞ przy $\mathbf {I} \;$ jako macierzy jednostkowej:

$||(\mathbf {I} +\mathbf {M} )^{-1}||_{p}\leq {{1} \over {1-||\mathbf {M} ||_{p}}}\;$

(5.15)

Dowód
Jeśli macierz $\mathbf {I} +\mathbf {M} \;$ byłaby macierzą osobliwą, to wtedy dla niezerowego $\mathbf {x} \;$ powinno być spełnione $(\mathbf {I} +\mathbf {M} )\mathbf {x} =0\;$, to wtedy z definicji normy powinien być spełniony warunek:

$||\mathbf {x} ||_{p}=||-\mathbf {M} \mathbf {x} ||_{p}\leq ||\mathbf {M} ||_{p}||\mathbf {x} ||_{p}\;$

(5.16)

stąd dla dowolnej macierzy $\mathbf {M} \;$ warunek (5.18) z warunkiem lematu $||\mathbf {M} ||_{p}<1\;$ nie jest spełniony, stąd otrzymaliśmy sprzeczność. Z definicji elementu odwrotnego możemy powiedzieć:

$(\mathbf {I} +\mathbf {M} )^{-1}(\mathbf {I} +\mathbf {M} )=\mathbf {I} \Rightarrow (\mathbf {I} +\mathbf {M} )^{-1}+(\mathbf {I} +\mathbf {M} )^{-1}\mathbf {M} =\mathbf {I} \;$

(5.17)

Z obliczeń napisanej w punkcie (5.17) możemy napisać następujący wniosek:

$1=||\mathbf {I} ||_{p}\geq ||(\mathbf {I} +\mathbf {M} )^{-1}||_{p}-||(\mathbf {I} +\mathbf {M} )^{-1}||_{p}||\mathbf {M} ||_{p}\Rightarrow ||(\mathbf {I} +\mathbf {M} )^{-1}||_{p}\leq {{1} \over {1-||\mathbf {M} ||_{p}}}\;$

(5.18)

Maksymalna wartość wartości własnej nazywamy promieniem spektralnym macierzy A i na podstawie wzoru (5.15) piszemy:

$\rho (A)\leq ||\mathbf {A} ||_{p}{\mbox{ dla }}p=1,2,\infty ,E\;$

(5.19)

Lemat drugi
Dla każdej macierzy i dla liczby ε istnieje indukowana norma $||\cdot ||_{p}\;$, dla której zachodzi:

$||\mathbf {A} ||_{p}\leq \rho (\mathbf {A} )+\epsilon \;$

(5.20)

Dowód
Niech λ będzie jedną z największej wartości własnej, wtedy z definicji norm podanej powyżej dla macierzy zdiagonalizowanej, ta norma jej jest równa wartości własnej największej rozważanej macierzy, co stąd jeśli $||\mathbf {A} ^{'}||_{pr}\;$ jest macierzą diagonalną, to możemy napisać z transformacji macierzy podobnych $\mathbf {A} ^{'}=\mathbf {U} ^{-1}\mathbf {A} \mathbf {U} \;$, to możemy powiedzieć:

$||\mathbf {U} ^{-1}||_{pr}|||\mathbf {A} ||_{pr}||U||_{pr}\geq ||A^{'}||\Rightarrow k||\mathbf {A} ||\geq ||\mathbf {A} ^{'}||$

(5.21)

Przeprowadźmy teraz dowód poprzez zaprzeczenie. Załózmy, że $||\mathbf {A} ||\leq ||\mathbf {A} ^{'}||\;$, wtedy da się wybrać takie k, by nierówność (5.21) była niespełniona, zatem $||\mathbf {A} ||\geq ||\mathbf {A} ^{'}||\;$. Analogicznie do (5.21) da się napisać tożsamość przy tym samym k:

$k||\mathbf {A} ^{'}||\geq ||\mathbf {A} ||\;$

(5.22)

Co przeprowadzając w dowodzie poprzez zaprzeczenie takie same wywody co dla (5.21) dochodzimy do wniosku $||\mathbf {A} ^{'}||\geq ||\mathbf {A} ||\;$, wtedy łącząc te dwie nierówności otrzymujemy, że normy macierzy podobnych są sobie równe.

Twierdzenie drugie
Ciąg wektorów $\mathbf {A} \mathbf {x} ,\mathbf {A} ^{2}\mathbf {x} ,...,\mathbf {A} ^{i}\mathbf {x} ,..$ jest zbieżny do zera, wtedy i tylko wtedy gdy promień spektralny macierzy jest mniejszy od jedynki.

Dowód
Udowodnijmy powyższe twierdzenie poprzez zaprzeczenie, zatem jeśli powyższy ciąg dąży zera, i promień spektralny jest większy lub równy niż zera, to powinno zachodzić na podstawie (5.21):

$0=\lim _{i\rightarrow \infty }||\mathbf {A} ^{i}\mathbf {x} ||\leq \lim _{i\rightarrow \infty }||\mathbf {A} ||_{p}^{i}||\mathbf {x} ||_{p}\geq \lim _{i\rightarrow \infty }\rho (\mathbf {A} )^{i}||\mathbf {x} ||_{p}\geq ||\mathbf {x} ||_{p}\;$

(5.23)

Ponieważ według (5.23) ciąg zależny od "i" dąży do zera i jednocześnie jest większy lub równy $||\mathbf {x} ||_{p}\;$ dla dowolnego $\mathbf {x} \;$, stad mamy sprzeczność. Udowodnijmy teraz twierdzenie odwrotne, jeśli ciąg jest rozbieżny od zera i promień spektralny jest mniejszy od zera dla dowolnie małego ε, to na podstawie wzoru (5.20) powiemy:

$\lim _{i\rightarrow \infty }||\mathbf {A} ^{i}\mathbf {x} ||\leq \lim _{i\rightarrow \infty }||\mathbf {A} ||_{p}^{i}||\mathbf {x} ||_{p}\leq \lim _{\overset {i\rightarrow \infty }{\underset {\epsilon \rightarrow 0}{}}}(\rho (\mathbf {A} )+\epsilon )^{i}=0\;$

(5.24)

Stąd jeśli rozważany ciąg jest zbieżny do zera według obliczeń (5.24), to nie jest jednocześnie rozbieżny, stąd sprzeczność, więc udowodniane twierdzenie jest prawdziwe.

Błędy rozwiązań układów równań algebraicznych

Przypadek pierwszy

Weźmy sobie zamiast $\mathbf {A} \;$ i $\mathbf {x} \;$ wielkości $\mathbf {A} +\delta \mathbf {A} \;$ i $\mathbf {x} +\delta \mathbf {x} \;$, wtedy równanie macierzowe $\mathbf {A} \mathbf {x} =\mathbf {b} \;$ możemy zapisać przy $\delta \mathbf {A} \;$ wzorem:

$\mathbf {A} (\mathbf {x} +\delta \mathbf {x} )=\mathbf {b} +\delta \mathbf {b} \Rightarrow \delta x=\mathbf {A} ^{-1}\delta \mathbf {b} \;$

(5.25)

Dla dowolnej normy wektora zaburzenia zapisany przy pomocy normy macierzy $\mathbf {A} \;$ dla równania macierzowego $\mathbf {A} \mathbf {x} =\mathbf {b} \;$ możemy zapisać $||\delta \mathbf {x} ||_{p}=||\mathbf {A} ^{-1}||_{pq}||\delta \mathbf {b} ||_{q}\;$, wtedy względne zaburzenie wielkości $\delta \mathbf {x} \;$ możemy zapisać przez względną zmianę wielkości $\mathbf {b} \;$, co daje nam w rezultacie przy oznaczeniu wielkości k:

${{||\delta \mathbf {x} ||_{p}} \over {||\mathbf {x} ||_{p}}}\leq \underbrace {{||\mathbf {A} ^{-1}||_{qp}||\mathbf {b} ||_{q}} \over {||\mathbf {x} ||_{p}}} _{k}{{||\delta \mathbf {b} ||_{q}} \over {||\mathbf {b} ||_{q}}}=k{{||\delta \mathbf {b} ||_{q}} \over {||\mathbf {b} ||_{q}}}\;$

(5.26)

wtedy w (5.26) zachodzi równość tylko dla ściśle określonego $\delta \mathbf {b} \;$, a nierówność słaba dla $\delta \mathbf {b} \;$, który należy od n-wymiarowej przestrzeni rzeczywistej. Jeśli skorzystamy z rozważanego równania macierzowego, wtedy definicję na czynnik "k" możemy przepisać w formie niezależącego od wektora zmiennych $\mathbf {x} \;$ w postaci pewnej nierówności:

$k={{||\mathbf {A} ^{-1}||_{pq}||\mathbf {A} \mathbf {x} ||_{q}} \over {||\mathbf {x} ||_{p}}}\leq ||\mathbf {A} ^{-1}||_{pq}||\mathbf {A} ||_{pq}=K_{pq}\;$

(5.27)

Liczbę K_pq będziemy nazywać wskaźnikiem uwarunkowania równania macierzowego $\mathbf {A} \mathbf {x} =\mathbf {b} \;$, który symbolizuje pewien układ równań.

Przypadek drugi

Niech dalej $\delta \mathbf {b} =0\;$ i $\delta \mathbf {A} \neq 0\;$, wtedy możemy napisać:

$(\mathbf {A} +\delta \mathbf {A} )(\mathbf {x} +\delta \mathbf {x} )=\mathbf {b} \Rightarrow \mathbf {A} \mathbf {x} +\mathbf {A} \delta \mathbf {x} +\delta \mathbf {A} \mathbf {x} +\delta \mathbf {A} \delta \mathbf {x} =\mathbf {b} \Rightarrow \mathbf {A} \delta \mathbf {x} +\delta \mathbf {A} \mathbf {x} =0\Rightarrow \;$
$\Rightarrow \delta \mathbf {x} =-\mathbf {A} ^{-1}\delta \mathbf {A} \mathbf {x} \Rightarrow ||\delta \mathbf {x} ||_{p}\leq ||\mathbf {A} ^{-1}||_{qp}||\delta \mathbf {A} ||_{pq}\cdot ||\mathbf {x} ||_{p}=\underbrace {||\mathbf {A} ^{-1}||_{qp}\cdot ||\mathbf {A} ||_{pq}} _{k}\cdot ||\mathbf {x} ||_{p}{{||\delta \mathbf {A} ||_{pq}} \over {||\mathbf {A} ||_{pq}}}\;$

(5.28)

Z obliczeń przeprowadzonych w punkcie (5.28) dochodzimy do wniosku przy odpowiedniej definicji k, że zachodzi:

${{||\delta \mathbf {x} ||_{p}} \over {||\mathbf {x} ||_{p}}}\leq k{{||\delta \mathbf {A} ||_{pq}} \over {||\mathbf {A} ||_{pq}}}\;$

(5.29)

Określmy jaka jest graniczna wartość współczynnika wzmocnienia k₀, jako wielkości granicznej współczynników k określonych we wzorze (5.28) jako czynnik, takiego że zachodzi równość (5.29) dla ściśle określonego $\delta \mathbf {A} \;$, która dla każdej tej wielkości zachodzi w ogólności nierówność słaba, takiego że:

${{||\delta \mathbf {A} ||} \over {||\mathbf {A} ||}}\leq \delta \;$

(5.30)

gdy zachodzi δ→0. Teraz będziemy wyznaczać oszacowanie wielkości k₀, które wyznaczymy poniżej Weźmy sobie do rozważenia równość, którą jest tożsamościowo równa zero, i ją zapisujemy sposobem: $F(\mathbf {x} ,\mathbf {A} )=\mathbf {A} \mathbf {x} -\mathbf {b} =0\;$. Weźmy sobie wektor $\mathbf {a} \;$, który przepiszemy w postaci: $\mathbf {a} =[a_{11},a_{21},...,a_{n1},a_{12},a_{22},...,a_{n2},..,a_{n1},a_{n2},...,a_{nn}]^{T}\;$. Mając sobie funkcję $F(\mathbf {x} ,\mathbf {a} )=0\;$, to wtedy na podstawie różniczki zupełnej tejże funkcji możemy powiedzieć:

$0={{\partial F(\mathbf {x} ,\mathbf {a} )} \over {\partial \mathbf {x} }}\delta \mathbf {x} +{{\partial F(\mathbf {x} ,\mathbf {a} )} \over {\partial \mathbf {a} }}\delta \mathbf {a} =F_{x}(\mathbf {x} ,\mathbf {a} )d\mathbf {a} +F_{a}(\mathbf {x} ,\mathbf {a} )\delta \mathbf {a} \Rightarrow {{\delta \mathbf {x} } \over {\delta \mathbf {a} }}=-[F_{x}(\mathbf {x} ,\mathbf {a} )]^{-1}F_{a}(\mathbf {x} ,\mathbf {a} )\;$

(5.31)

wtedy możemy napisać wielkość $F_{x}(\mathbf {x} ,\mathbf {a} )=\mathbf {A} \;$, która jest macierzą nieosobliwą, a $F_{x}(\mathbf {x} ,\mathbf {a} )=[x_{1}\mathbf {I} ,x_{2}\mathbf {I} ,...,x_{n}\mathbf {I} ]\;$, zatem po wykorzystaniu tożsamości na $F_{a}(\mathbf {x} ,\mathbf {a} )\;$ i na $F_{x}(\mathbf {x} ,\mathbf {a} )\;$, wtedy różniczkę wielkości $\mathbf {x} \;$ możemy przepisać w następującej formie:

$\delta \mathbf {x} =-[x_{1}\mathbf {A} ^{-1},x_{2}\mathbf {A} ^{-1},...,x_{n}\mathbf {A} ^{-1}]\delta \mathbf {a} +O(||\delta \mathbf {a} ||)\;$

(5.32)

Z definicji normy możemy napisać oszacowanie czynnika stojącego przy $\delta \mathbf {a} \;$ występującego we wzorze (5.32):

$||[x_{1}\mathbf {A} ^{-1},x_{2}\mathbf {A} ^{-1},..,x_{n}\mathbf {A} ^{-1}]||_{1\infty }=\max _{i,j,k}|x_{k}a_{ijk}^{-1}|=||\mathbf {x} ||_{\infty }||\mathbf {A} ^{-1}||_{1\infty }\;$

(5.33)

stąd możemy napisać oszacowanie wielkości ${{||\delta \mathbf {x} ||_{\infty }} \over {||\mathbf {x} ||_{\infty }}}\;$ poprzez wielkość ${{||\delta \mathbf {a} ||_{\infty }} \over {||\mathbf {A} ||_{1\infty }}}\;$, zatem:

${{||\delta \mathbf {x} ||_{\infty }} \over {||\mathbf {x} ||_{\infty }}}\leq \underbrace {||\mathbf {A} ^{-1}||_{1\infty }||\mathbf {A} ||_{1\infty }} _{k_{0}}{{||\delta \mathbf {a} ||_{1}} \over {||\mathbf {A} ||_{1\infty }}}+O(||\delta \mathbf {a} ||)\;$

(5.34)

Czynnik stojący przy k₀ możemy wyrazić przy pomocy definicji normy wektora (5.1), a także normy macierzy (5.8) w sposób ${{\delta \mathbf {a} ||_{1}} \over {||\mathbf {A} ||_{1\infty }}}={{\sum _{i=1}^{n}\sum _{j=1}^{n}|\delta a_{ij}|} \over {\max _{ij}|a_{ij}|}}\;$, a k₀ oblicza się jako iloczyn norm macierzy macierzy A i jego odwrotności przy normie (5.8), który jest iloczynem maksymalnych wartości modułów elementów macierzowych macierzy A i jego odwrotności.

Przypadek trzeci

Weźmy sobie, że różniczki wyrazów wolnych $\mathbf {b} \;$ i macierzy $\mathbf {A} \;$, które są nierówne zero, a macierz $\mathbf {A} +\delta \mathbf {A} \;$ jest macierzą nieosobliwą, wtedy:

$(\mathbf {A} +\delta \mathbf {A} )(\mathbf {x} +\delta \mathbf {x} )=\mathbf {b} +\delta \mathbf {b} \Rightarrow (\mathbf {A} +\delta \mathbf {A} )\mathbf {x} +(\mathbf {A} +\delta \mathbf {A} )\delta \mathbf {x} =\mathbf {b} +\delta \mathbf {b} \Rightarrow$
$\Rightarrow \delta \mathbf {A} \mathbf {x} +(\mathbf {A} +\delta \mathbf {A} )\delta \mathbf {x} =\delta \mathbf {b} \Rightarrow \delta \mathbf {x} =\left(1+\mathbf {A} ^{-1}\delta \mathbf {A} \right)^{-1}\mathbf {A} ^{-1}\left(\delta \mathbf {b} -\delta \mathbf {A} \mathbf {x} \right)\;$

(5.35)

Względny błąd wielkości $\delta \mathbf {x} \;$ możemy napisać po wykorzystaniu równania macierzowego $\mathbf {A} \mathbf {x} =\mathbf {b} \;$, którą wykorzystamy do wyliczenia nierówności drugiej poniżej występującego w nawiasie po prawej stronie jako drugi czynnik w pierwszym składniku, które jest zawsze mniejsze od jedynki według dowodu ${{||\mathbf {b} ||_{q}} \over {||\mathbf {x} ||_{p}||\mathbf {A} ||_{pq}}}={{||\mathbf {A} \mathbf {x} ||_{q}} \over {||\mathbf {x} ||_{p}||\mathbf {A} ||_{pq}}}\leq 1\;$, a także z małości zaburzenia $\delta \mathbf {A} \;$ będziemy mogli zapisać $||\mathbf {A} ^{-1}||\cdot ||\delta \mathbf {A} ||=\rho <1\;$, w takim razie:

${{||\delta \mathbf {x} ||_{p}} \over {||\mathbf {x} ||_{p}}}\leq {{1} \over {1-||\mathbf {A} ^{-1}\delta \mathbf {A} ||_{pq}}}||\mathbf {A} ^{-1}||_{qp}\cdot ||\mathbf {A} ||_{pq}\left({{||\delta \mathbf {b} ||_{q}} \over {||\mathbf {A} ||_{pq}\cdot ||\mathbf {x} ||_{p}}}+{{||\delta \mathbf {A} ||_{pq}} \over {||\mathbf {A} ||_{pq}}}\right)\Rightarrow \;$
$\Rightarrow {{||\delta \mathbf {x} ||_{p}} \over {||\mathbf {x} ||_{p}}}\leq {{1} \over {1-||\mathbf {A} ^{-1}||_{qp}\cdot ||\delta \mathbf {A} ||_{pq}}}||\mathbf {A} ^{-1}||_{qp}\cdot ||\mathbf {A} ||_{pq}\left({{||\delta \mathbf {b} ||_{q}} \over {||\mathbf {b} ||_{q}}}{{||\mathbf {b} ||_{q}} \over {||\mathbf {x} ||_{p}||\mathbf {A} ||_{pq}}}+{{||\delta \mathbf {A} ||_{pq}} \over {||\mathbf {A} ||_{pq}}}\right)\Rightarrow \;$
$\Rightarrow {{||\delta \mathbf {x} ||_{p}} \over {||\mathbf {x} ||_{p}}}\leq {{1} \over {1-||\mathbf {A} ^{-1}||_{qp}\cdot ||\delta \mathbf {A} ||_{pq}}}||\mathbf {A} ^{-1}||_{qp}\cdot ||\mathbf {A} ||_{pq}\left({{||\delta \mathbf {b} ||_{q}} \over {||\mathbf {b} ||_{q}}}+{{||\delta \mathbf {A} ||_{pq}} \over {||\mathbf {A} ||_{pq}}}\right)\Rightarrow \;$
$\Rightarrow {{||\delta \mathbf {x} ||_{p}} \over {||\mathbf {x} ||_{p}}}\leq {{1} \over {1-\rho }}||\mathbf {A} ^{-1}||_{qp}\cdot ||\mathbf {A} ||_{pq}\left({{||\delta \mathbf {b} ||_{q}} \over {||\mathbf {b} ||_{q}}}+{{||\delta \mathbf {A} ||_{pq}} \over {||\mathbf {A} ||_{pq}}}\right)\;$

(5.36)

Układy równań algebraicznych o trójkątnej macierzy

Weźmy sobie układ równań, w której wszystkie wyrazy na diagonalnej są różne od zera, to możemy wyznaczyć zmienną $\mathbf {x} \;$ wykonując małą liczbę działań arytmetycznych i dla małych błędów zaokrągleń przy wyznaczaniu poszczególnych elementów wspomnianego wektora, zatem ten nasz układ równań piszemy:

a₁₁x₁+a₁₂x₂+...+a_1nx_n=b₁

..........a₂₂x₂+...+a_2nx_n=b₂

........................................

.........................a_nnx_n=b_n

(5.37)

Poszczególne elementy zmiennych x_k możemy policzyć z poniższych wzorów:

$x_{n}={{b_{n}} \over {a_{nn}}}\;$

$x_{n-1}={{b_{n-1}-a_{n-1n}x_{n}} \over {a_{n-1n-1}}}\;$

$x_{i}={{b_{i}-a_{in}x_{n}-...-a_{ii+1}x_{i+1}} \over {a_{ii}}}{\mbox{ dla }}i=n-1,n-2,..,1\;$

(5.38)

Patrząc na powyższy układ równań można powiedzieć, że dla i=n jest wykonywane jedno dzielenie lub mnożenie, a dla i=n-1 jest wykonywane jedno mnożenie i jedno dzielenie, czyli dwa mnożenia i dzielenia, a przy x₁ jest wykonywanych n dzieleń i mnożeń, wtedy ilość dzieleń i mnożeń jest:

$M=1+2+3+..n={{1+n} \over {2}}n={{1} \over {2}}n^{2}+{{1} \over {2}}n\;$

(5.39)

Dla układu rozwiązań x_i (5.34) dla x_n mamy zero dodawań i odejmowań, dla x_n-1 mamy jedno dodawanie lub odejmowanie, dla x₁ mamy n-1 dodawań lub odejmowań, zatem liczba wszystkich dodawań lub odejmowań jakie należy wykonać jest:

$D=0+1+2+...+n-1={{0+n-1} \over {2}}n={{1} \over {2}}n^{2}-{{1} \over {2}}n\;$

(5.40)

Jeśli będziemy wykonywali działania na komputerze, jeśli nie możemy wsadzić w miejsce wektora $\mathbf {b} \;$ wektora $\mathbf {x} \;$, to liczba komórek zajętych dla i=n jest 3, dla i=n-1 jest 4, i idać dalej dochodzimy aż do i=1 jest n+2 zajętych komórek pamięci, zatem liczba wszystkich komórek zajętych jest:

$P=3+4+...n+2={{3+n+2} \over {2}}n={{n+5} \over {2}}n={{1} \over {2}}n^{2}+{{5} \over {2}}n\;$

(5.41)

Wzory (5.38) możemy napisać z odpowiednim przybliżeniem uwzględniając, że wielkości x_i i b_i są napisane z pewnym zaokrągleniem:

$|\delta a_{ij}|\leq \epsilon {\begin{bmatrix}n|a_{11}|&(n+2)|a_{1,2}|&\cdots &4|a_{1,n-1}|&3|a_{1.n}|\\0&(n-1)|a_{2,2}|&\cdots &4|a_{1,n-1}|&3|a_{1.n}|\\\vdots &\vdots &\ddots &\vdots &\vdots \\0&0&\cdots &2|a_{n-1,n-1}&3|a_{n-1,n}|\\0&0&\cdots &0&1|a_{nn}|\end{bmatrix}}\;$

(5.42)

Powyższa macierzowa "nierówność" stanowi układ 1+2+3+..+n=n(n+1)/2 nierówności. Z powyższego oszacowania możemy wyedukować, że dla dowolnej normy p=1,∞,E spełniona jest zawsze nierówność:

$||\delta \mathbf {A} ||_{p}\leq \epsilon (n+2)||\mathbf {A} ||_{p}\;$

(5.43)

Jeśli wiadomo, że oszacowanie elementów macierzy $||\mathbf {A} ||\;$ z definicji normy ||⋅||_1∞ od góry jest $||\mathbf {A} ||_{1\infty }\leq g\;$, wtedy macierz błędu z normą ||⋅||₁ i ||⋅||_∞ możemy przepisać:

$||\delta \mathbf {A} ||_{1}\leq \epsilon \left({{1} \over {4}}n^{2}+{{5} \over {2}}n\right)g\;$

(5.44)

$||\delta \mathbf {A} ||_{\infty }\leq \epsilon \left({{1} \over {2}}n^{2}+{{5} \over {2}}n-2\right)g\;$

(5.45)

Wprowadźmy definicję macierzy A diagonalnie dominującej, których moduły elementów na diagonali są większe od sumy elementów macierzy stojącej w danym wierszu bez elementu na diagonali. Jest to macierz spełniającej warunek dla i=1,2,...,n:

$|a_{ii}|\geq \sum _{\overset {k=1}{\underset {k\neq i}{}}}^{n}|a_{ik}|\;$

(5.46)

Macierzą silnie diagonalnie dominującą nazywamy taką macierz A o stopniu n, w której w (5.46) występuje nierówność ostra. Macierz (silnie) diagonalnie dominująca kolumnowo nazywamy macierz, dla której A^T jest (silnie) diagonalnie dominująca, tzn. gdy spełniony jest warunek dla i=1,2,...,n:

$|a_{ii}|\geq (>)\sum _{\overset {k=1}{\underset {k\neq i}{}}}^{n}|a_{ki}|\;$

(5.47)

Macierz A generowana przez (5.37) nazywamy diagonalnie dominującą, gdy spełniony jest warunek pierwszy poniżej, a także ta macierz jest diagonalnie dominująca kolumnowo, gdy spełniony jest drugi warunek poniżej:

$||\delta \mathbf {A} ^{(z)}||_{\infty }\leq \epsilon (2n+2)g\;$

(5.48)

$||\delta \mathbf {A} ^{(z)}||_{1}\leq \epsilon (2n+1)g\;$

(5.49)

Rozwiązania równań liniowych metodą eliminacji Gaussa

Będziemy rozwiązywać układ równań liniowych z pełną macierzą $\mathbf {A} \;$, i doprowadzać ten układ równań do postaci trójkątnej, czyli przy wykorzystaniu równań (5.38), by później obliczyć jej poszczególne zmienne. Na sam początek określmy równanie macierzowe $\mathbf {A} ^{(1)}\mathbf {b} =\mathbf {b} ^{(1)}\;$, które możemy przepisać w postaci układu równań:

$a_{11}^{(1)}x_{1}+a_{12}^{(1)}x_{2}+...+a_{1n}^{(1)}x_{n}=b_{1}^{(1)}\;$

$a_{21}^{(1)}x_{1}+a_{22}^{(1)}x_{2}+...+a_{2n}^{(1)}x_{n}=b_{2}^{(1)}\;$
..........................................................

$a_{n1}^{(1)}x_{1}+a_{n2}^{(1)}x_{2}+...+a_{nn}^{(1)}x_{n}=b_{n}^{(1)}\;$

(5.50)

Odejmijmy i-te równanie i>1 od pierwszego równania pomnożonego przez $a_{i1}^{(1)}/a_{11}^{(1)}\;$ układu równań (5.46), w ten sposób otrzymujemy równanie macierzowe $\mathbf {A} ^{(2)}\mathbf {x} =\mathbf {b} ^{(2)}\;$, co po rozpisaniu jego na układ równań, otrzymujemy:

$a_{11}^{(2)}x_{1}+a_{12}^{(2)}x_{2}+...+a_{1n}^{(2)}x_{n}=b_{1}^{(2)}\;$

$a_{22}^{(2)}x_{2}+...+a_{2n}^{(2)}x_{n}=b_{2}^{(2)}\;$
..........................................................

$a_{n2}^{(2)}x_{2}+...+a_{nn}^{(2)}x_{n}=b_{n}^{(2)}\;$

(5.51)

W powyższych równaniu wyeliminowaliśmy zmienną x₁ dla równania o numerze więcej niż pierwszy. Idąc dalej tym sposobem od równania i>2 odejmujemy od równania drugiego pomnożonej przez $a_{i2}^{(2)}/a_{22}^{(1)}\;$, w ten sposób otrzymujemy układ macierzowy $\mathbf {A} ^{(3)}\mathbf {x} =\mathbf {b} ^{(2)}\;$. Dokonując tą metodę dalej według schematu podanego powyżej otrzymujemy na samym końcu równanie macierzowe $\mathbf {A} ^{(n)}\mathbf {x} =\mathbf {b} ^{(n)}\;$, które zapisujemy w postaci układu równań:

$a_{11}^{(n)}x_{1}+a_{12}^{(n)}x_{2}+...+a_{1n}^{(2)}x_{n}=b_{1}^{(n)}\;$

$a_{22}^{(n)}x_{2}+...+a_{2n}^{(n)}x_{n}=b_{2}^{(n)}\;$
..........................................................

$a_{nn}^{(n)}x_{n}=b_{n}^{(n)}\;$

(5.52)

W uzyskanym układzie równań (5.52) rozwiązujemy x_i metodą (5.38) dla macierzy $\mathbf {A} \;$ trójkątnej. Liczba mnożeń i dzieleń jakie należy wykonać dla układu równań (5.50), by uzyskać macierz trójkątna $\mathbf {A} \;$ dla układu równań (5.48) jest w ilości kroków:

$D=\sum _{i=0}^{n-1}(n+1-i)(n-1-i)=\sum _{i=1}^{n}(i+1)(i-1)=\sum _{i=1}^{n}i^{2}-\sum _{i=1}^{n}1={{n^{3}} \over {3}}+{{n^{2}} \over {2}}+{{n} \over {6}}-n=\;$
$={{n^{3}} \over {3}}+{{1} \over {2}}n^{2}-{{5} \over {6}}n\;$

(5.53)

Powyższy rozwiązanie należy uzyskać, jeśli dla pierwszego kroku należy wykonać dzielenie $a_{i1}^{(1)}/a_{11}^{(1)}\;$, które jest jedno, i wymnożenie jej przez n wyrazów pierwszego równania układu równań (5.50) poczynając od drugiego (bo pierwszy składnik chociaż jest wymnażany przez tą liczbę, ale tylko służy do wyzerowania pierwszego składnika w równania drugiego), tych dzieleń i mnożeń w sumie jest n+1, dla i=2,..,n, to aby przejść do układu równań (5.51) należy wykonać dodawań i mnożeń w ilości (n+1)(n-1), aby przejść od układu równań oznaczonej "i" na do i+1 należy wykonać działań (n+1-i)(n-1-i), czyli w sumie mamy takich naszych działań (5.53). Liczba dodawań i dzieleń jakie układ musi wykonać by układ (5.50) doprowadzić do postaci trójkątnej jest:

$D=\sum _{i=0}^{n-1}(n-i)(n-1-i)=\sum _{i=1}^{n}i(i-1)=\sum _{i=1}^{n}i^{2}-\sum _{i=1}^{n}i={{n^{3}} \over {3}}+{{n^{2}} \over {2}}+{{n} \over {6}}-{{n} \over {2}}-{{n^{2}} \over {2}}=\;$
$={{n^{3}} \over {3}}-{{1} \over {3}}n\;$

(5.54)

Powyższe rozwiązanie możemy uzyskać, gdy dla pierwszego układu równań przejdziemy do drugiego układu równań, to wtedy należy wykonać n(n-1) działań, to aby przejść z układu równań i-tego do i+1-ego, to należy wykonać (n-i)(n-1-i) dodawań, odejmowań, zatem całkowita ilość dodawań i odejmowań aby układ był w postaci trójkątnej jest uzyskana równaniem (5.54). Powyższy sposób przeprowadziliśmy bez ujawniania jawnego operacji macierzowych, tylko powiedzieliśmy co dokonaliśmy na poszczególnych równaniach, do tego problemu podejdźmy z innej strony, tzn. za pomocą operacji macierzowych, aby przejść do równoważnych równań macierzowych z $\mathbf {A} ^{(1)}\mathbf {x} =\mathbf {b} ^{(1)}\;$ do $\mathbf {A} ^{(2)}\mathbf {x} =\mathbf {b} ^{(2)}\;$, to pierwsze z tych równań macierzowych należy pomnożyć obustronnie lewostronnie przez macierz:

$\mathbf {L} ^{(1)}={\begin{bmatrix}1&0&0&\vdots &0\\-l_{21}&0&0&\vdots &0\\\vdots &\vdots &\vdots &\ddots &\vdots \\-l_{n1}&0&0&\cdots &1\\\end{bmatrix}}{\mbox{ dla }}l_{i1}={{a_{i1}^{(1)}} \over {a_{11}^{(1)}}}{\mbox{, }}i=2,3,..n\;$

(5.55)

Przy pomocy macierzy podanej w podpunkcie (5.55) możemy pomnożyć początkowe równanie macierzowe lewostronnie $\mathbf {A} ^{(1)}\mathbf {x} =\mathbf {b} ^{(1)}\;$, w ten sposób uzyskując równanie macierzowe $\mathbf {L} ^{(1)}\mathbf {A} ^{(1)}\mathbf {x} =\mathbf {L} ^{(1)}\mathbf {b} ^{(1)}\;$. By otrzymać równość macierzową $\mathbf {L} ^{(3)}\mathbf {x} =\mathbf {b} ^{(3)}\;$, to drugie równanie macierzowe należy pomnożyć obustronnie lewostronnie przez macierz:

$\mathbf {L} ^{(2)}={\begin{bmatrix}1&0&0&\cdots &0\\0&1&0&\cdots &0\\0&-l_{32}&1&\cdots &0\\\vdots &\vdots &\vdots &\ddots &\vdots \\0&-l_{n2}&0&\cdots &1\end{bmatrix}}{\mbox{ dla }}l_{i1}={{a_{i2}^{(1)}} \over {a_{22}^{(1)}}}{\mbox{, }}i=3,4,..n\;$

(5.56)

Macierz końcową $\mathbf {A} ^{(n)}\;$, która jest macierzą trójkątną, to by go otrzymać należy macierz $\mathbf {A} ^{(1)}\;$ wymnożyć lewostronnie przez macierze $\mathbf {L} ^{(i)}\;$, gdzie i=2,3,..n, a także to samo robimy uzyskując wektor wyrazów wolnych $\mathbf {b} ^{(n)}\;$, co je napiszemy w jednej linijce:

$\mathbf {A} ^{(n)}=\mathbf {L} ^{(n-1)}\mathbf {L} ^{(n-2)}...\mathbf {L} ^{(1)}\mathbf {A} ^{(1)}\;$

(5.57)

$\mathbf {b} ^{(n)}=\mathbf {L} ^{(n-1)}\mathbf {L} ^{(n-2)}...\mathbf {L} ^{(1)}\mathbf {b} ^{(1)}\;$

(5.58)

Ponieważ macierze $\mathbf {L} ^{i)}\;$ są macierzami nieosobliwymi, wtedy możemy równanie (5.57) równoważnie zapisać na w sposób:

$\mathbf {A} ^{(1)}=(\mathbf {L} ^{(1)})^{-1}(\mathbf {L} ^{(2)})^{-1}...(\mathbf {L} ^{(n-1)})^{-1}\mathbf {A} ^{(n)}\;$

(5.59)

Można udowodnić, że odwrotności macierzy (5.55) i (5.56), czyli macierze $(\mathbf {L} ^{(1)})^{-1}\;$ i $(\mathbf {L} ^{(2)})^{-2}\;$, zapisujemy w formie:

$(\mathbf {L} ^{(1)})^{-1}={\begin{bmatrix}1&0&0&\vdots &0\\l_{21}&0&0&\vdots &0\\\vdots &\vdots &\vdots &\ddots &\vdots \\l_{n1}&0&0&\cdots &1\\\end{bmatrix}}\;$

(5.60)

$(\mathbf {L} ^{(2)})^{-1}={\begin{bmatrix}1&0&0&\cdots &0\\0&1&0&\cdots &0\\0&l_{32}&1&\cdots &0\\\vdots &\vdots &\vdots &\ddots &\vdots \\0&l_{n2}&0&\cdots &1\end{bmatrix}}\;$

(5.61)

Możemy również powiedzieć, że iloczyn odwrotności macierzy $\mathbf {L} ^{(i)}\;$ dla i=1,...,n-1 możemy zapisać przy pomocy liczb l_ij w formie:

$(\mathbf {L} ^{(1)})^{-1}(\mathbf {L} ^{(2)})^{-1}...(\mathbf {L} ^{(n-1)})^{-1}={\begin{bmatrix}1&0&0&\cdots &0\\l_{21}&1&0&\cdots &0\\l_{31}&l_{32}&1&\cdots &0\\\vdots &\vdots &\vdots &\ddots &\vdots \\l_{n1}&l_{n2}&l_{n3}&\cdots &1\end{bmatrix}}\;$

(5.62)

Oznaczmy macierz oznaczoną (5.62) przez $\mathbf {L} \;$, a macierz $\mathbf {A} ^{(n)}\;$ przez $\mathbf {U} \;$, wtedy wyrażenie macierzowe (5.59) zapisujemy jako:

$\mathbf {A} ^{(1)}=\mathbf {L} \mathbf {U} \;$

(5.63)

Macierz $\mathbf {L} \;$ przedstawiona powyżej jest macierzą trójkątną dolną. Równanie macierzowe $\mathbf {A} \mathbf {x} =\mathbf {b} \;$ zapisujemy w formie $\mathbf {L} \mathbf {U} \mathbf {x} =\mathbf {b} \;$, co można zapisać po przekształceniu $\mathbf {U} \mathbf {x} =\mathbf {L} ^{-1}\mathbf {b} =\mathbf {b} ^{(n)}\;$, co jest równoważne rozwiązaniu układu równań macierzowych $\mathbf {L} \mathbf {y} =\mathbf {b} \;$ i $\mathbf {U} \mathbf {x} =\mathbf {y} \;$, z których będziemy szukali rozwiązania $\mathbf {x} \;$. Jeśli znamy macierz LU, to ilość mnożeń jakich należy dokonać przy wyznaczeniu wektora "x" jest wyrażona przez $M={{1} \over {2}}n(n-1)+{{1} \over {2}}n(n+1)=n^{2}\;$, a ilość dodawań jakie musimy dokonać jest napisana przez $D=n^{2}-n\;$.

Wyznaczanie macierzy L i U, a metoda Doolittle'a

Równanie macierzowe (5.63) możemy przepisać w postaci rozwiniętej w takiej postaci, tzn. w której napiszemy zamiast symboli macierzy wchodzących w skład tego równania ich postać rozwiniętą obrazująca poszczególne ich elementy:

${\begin{bmatrix}a_{11}&a_{12}&\cdots &a_{1n}\\a_{21}&a_{22}&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &a_{nn}\end{bmatrix}}={\begin{bmatrix}1&0&\cdots &0\\l_{21}&1&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\l_{n1}&l_{n2}&\cdots &1\end{bmatrix}}\cdot {\begin{bmatrix}u_{11}&u_{12}&\cdots &u_{1n}\\0&u_{22}&\cdots &u_{2n}\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &u_{nn}\end{bmatrix}}\;$

(5.64)

Jeśli popatrzymy na równość macierzową (5.64) i wykonamy działanie w lewej jego stronie, to otrzymamy równość skalarną na elementy a_ij w których chcemy wyznaczyć l_ij i u_ij:

$a_{ij}=u_{ij}+\sum _{k=1}^{i-1}l_{ik}u_{kj}\Rightarrow u_{ij}=a_{ij}-\sum _{k=1}^{i-1}l_{ik}u_{kj}{\mbox{, dla }}j=i,i+1,...,n\;$

(5.65)

Równość (5.64) możemy przetransponować obustronnie, tzn. zamieniając wiersze z kolumnami, co w końcu otrzymujemy równość macierzową:

${\begin{bmatrix}a_{11}&a_{21}&\cdots &a_{n1}\\a_{12}&a_{22}&\cdots &a_{n2}\\\vdots &\vdots &\ddots &\vdots \\a_{1n}&a_{2n}&\cdots &a_{nn}\end{bmatrix}}={\begin{bmatrix}u_{11}&0&\cdots &0\\u_{12}&u_{22}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\u_{1n}&u_{2n}&\cdots &u_{nn}\end{bmatrix}}\cdot {\begin{bmatrix}1&l_{12}&\cdots &l_{1n}\\0&1&\cdots &l_{2n}\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &1\end{bmatrix}}\;$

(5.66)

Mając do dyspozycji równość macierzową (5.66) możemy go przepisać ujmując poszczególne jej elementy dla j=i+1,i+2,...,n:

$a_{ji}=\sum _{k=1}^{i}u_{ki}l_{jk}\Rightarrow a_{ji}=u_{ii}l_{ji}+\sum _{k=1}^{i-1}l_{jk}u_{ki}\Rightarrow l_{ji}=\left(a_{ji}-\sum _{k=1}^{i-1}l_{jk}u_{ki}\right)/u_{ii}\;$

(5.67)

Ilość mnożeń w równaniach (5.65) i (5.67) jest wyrażona $M={{1} \over {3}}n^{3}-{{1} \over {3}}n\;$, a liczba dodawań w tych samych równaniach co poprzednio jest $D={{1} \over {3}}n^{3}-{{1} \over {2}}n^{2}+{{1} \over {6}}n\;$. Znając ile należy dokonać działań aby wyznaczyć wektor "x", gdy znamy macierz LU, to dodając te ilości do poprzednich M i D otrzymujemy liczbę dodawań i mnożeń taką samą jak przy metodzie Gaussa, tzn. ostateczną ilość mnożeń jaką musimy dokonać przy wyznaczeniu macierzy L i U, a potem do wyznaczania wektora "x" jest napisana $M={{1} \over {3}}n^{3}+n^{2}-{{1} \over {3}}n\;$, a liczba dodawań $D={{1} \over {3}}n^{3}+{{1} \over {2}}n^{2}-{{5} \over {6}}n\;$.

Zastosujmy zmodyfikowaną metodę eliminacją zwaną również częściowym wyborem elementu podstawowego. W tej metodzie elementem podstawowym nazywamy taki element macierzy A przy pomocy której eliminujemy zmienną w pozostałych równaniach. Te elementy zwykle przyjmuje się jako elementy diagonalne macierzy podstawowych A^(k), gdzie k=1,2,..,n. Te elementy wybieramy w taki sposób by one leżały w k-tek kolumnie w k-tego macierzy, by ten elementem miał największy moduł w stosunku do pozostałych elementów w danym wierszu macierzy A, co można je uzyskać poprzestawiając poszczególne wiersze w macierzy A, by później leżały one na diagonali, oczywiste jest, że również musimy jednocześnie przedstawiać wiersze wektorów x i b. Ta wersja algorytmu Gaussa z częściowym wyborem elementy częściowego nazywamy metodą Gaussa-Crouta, którą w dalszym ciągu będziemy nazywać metodą GCW. Ta metoda gwarantuje, że proces poszukiwania x nie zatrzyma się na dzieleniu przez zero przy założeniu nieosobliwości macierzy A.

Błędy przybliżeń w metodzie Gaussa i Doolittle'a

Obie metody tutaj rozważane są ze sobą równoważne numeryczne o tej samej precyzji, aby wyznaczyć wektor "x" należy wyznaczyć macierz ${\tilde {A}}=\mathbf {L} \mathbf {U} \;$ metodą Gaussa lub metodą Doolittle'a, a następnie znaleźć wektor "y" w równaniu macierzowym $\mathbf {L} \mathbf {y} =\mathbf {b} \;$, i dalej wyznaczyć "x" w równaniu $\mathbf {U} \mathbf {x} =\mathbf {y} \;$. Ponieważ macierze L i U są napisane z pewnym zaokrągleniem, to błąd zaokrąglenia LU napiszmy przez E, która jest błędem rozkładu i wtedy możemy powiedzieć $\mathbf {L} \mathbf {U} ={\tilde {\mathbf {A} }}+\mathbf {E} \;$. Znając błędu rozkładu L i U możemy napisać równości $(\mathbf {L} +\delta \mathbf {L} )\mathbf {x} =\mathbf {b} \;$ i $(\mathbf {U} +\delta \mathbf {U} )\mathbf {x} =\mathbf {y} \;$, co na podstawie tego możemy powiedzieć

$(\mathbf {L} +\delta \mathbf {L} )(\mathbf {U} +\delta \mathbf {U} )\mathbf {x} =\mathbf {b} \Rightarrow \left(\mathbf {L} \mathbf {U} +\mathbf {L} \delta \mathbf {U} +\delta \mathbf {L} \mathbf {U} +\delta \mathbf {L} \delta \mathbf {U} \right)\mathbf {x} =\mathbf {b} \;$

(5.68)

Z równości (5.68) możemy wyznaczyć błąd zaokrągleń macierzy A wiedząc jakie jest błąd zaokrąglenia macierzy E, zatem:

$\delta \mathbf {A} ^{(z)}=\mathbf {E} +\mathbf {L} \delta \mathbf {U} +\delta \mathbf {L} \mathbf {U} +\delta \mathbf {L} \delta \mathbf {U} \;$

(5.69)

Jeśli będziemy przyjmować $g_{ij}=\max _{k=1,2,..,n}|a_{ij}^{(k)}|\;$, to oszacowanie błędu macierzy LU dla poszczególnych jego elementów jest |e_ij|≤3,01εng_ij. Normy błędu zaokrągleń macierzy LU w metodzie GCW, pamiętając przy okazji, że g jest to największy moduł elementów macierzy ${\tilde {\mathbf {A} }}=\mathbf {A} ^{(1)},\mathbf {A} ^{(2)},...,\mathbf {A} ^{(n)}\;$, określamy według:

$||\mathbf {E} ||_{1}\leq \epsilon (n^{2}-n)g\;$

(5.70)

$||\mathbf {E} ||_{\infty }\leq \epsilon (n^{2}-n-2)g\;$

(5.71)

W przypadku metody GCW wielkość g spełnia nierówność g≤2^n-1||A||_1∞, ale jak się okazuje, że dobrym ograniczeniem dla g od góry jest nierówność $g\leq 8||{\tilde {\mathbf {A} }}||_{1\infty }\;$, a ponadto zachodzą również warunki $||{\tilde {\mathbf {A} }}||_{1}\geq ||{\tilde {\mathbf {A} }}||_{1\infty }\;$, $||{\tilde {\mathbf {A} }}||_{\infty }\geq ||{\tilde {\mathbf {A} }}||_{1\infty }\;$, wtedy na podstawie (5.70) i (5.71) można przejąć względne błędy macierzy A według oszacowań:

${{||\mathbf {E} ||_{1}} \over {||\mathbf {A} ||_{1}}}\leq 8\epsilon (n^{2}-n)\;$

(5.72)

${{||\mathbf {E} ||_{\infty }} \over {||\mathbf {A} ||_{\infty }}}\leq 8\epsilon (n^{2}+n-2)\;$

(5.73)

Jeśli dodatkowo uwzględnimy, że moduł poszczególnych elementów macierzy L jest mniejszy lub równy jeden, wtedy normy błędów macierzy L przepisujemy w formie:

$||\delta \mathbf {L} ||_{1}\leq \epsilon \left({{1} \over {4}}n^{2}+{{5} \over {2}}n-1\right)\;$

(5.74)

$||\delta \mathbf {L} ||_{\infty }\leq \epsilon \left({{1} \over {2}}n^{2}+{{5} \over {2}}n-3\right)\;$

(5.75)

A ilorazy błędu U przez jakąś normę A możemy przepisać w formie zależnego od stopnia n macierzy i zmiennej dowolnej ε:

${{||\delta \mathbf {U} ||_{1}} \over {||\mathbf {A} ||_{1}}}\leq \left({{1} \over {4}}n^{2}+{{5} \over {2}}n\right)\;$

(5.76)

${{||\delta \mathbf {U} ||_{\infty }} \over {||\mathbf {A} ||_{\infty }}}\leq \epsilon \left({{1} \over {2}}n^{2}+{{5} \over {2}}n-2\right)\;$

(5.77)

Z definicji normy ||⋅||₁ i ||⋅||_∞ możemy napisać oszacowania dla norm macierzy U i L przy definicji parametru g i stopnia wspomnianych macierzy, które zapisujemy w formie ||U||₁≤ng, ||U||_∞≤ng, ||L||_∞≤n, wtedy norma błędu bezwzględnego zaokrąglenia macierzy A według powyższych uwag przedstawiamy jako:

${{||\delta \mathbf {A} ^{(z)}||_{\infty }} \over {||\mathbf {A} ||_{\infty }}}\leq {{||\mathbf {E} ||_{\infty }} \over {||{\tilde {\mathbf {A} }}||_{\infty }}}+||\mathbf {L} ||_{\infty }{{||\delta \mathbf {U} ||_{\infty }} \over {||{\tilde {\mathbf {A} }}||_{\infty }}}+||\delta \mathbf {L} ||_{\infty }{{||\mathbf {U} ||_{\infty }} \over {||{\tilde {\mathbf {A} }}||_{\infty }}}+o(\epsilon )=\;$
$=\epsilon \left({{9} \over {2}}n^{3}+{{61} \over {2}}n^{2}-18n-16\right)+o(\epsilon )$

(5.78)

Jeśli popatrzymy na wzór (5.36) i oznaczymy $K_{\infty }=||\mathbf {A} ||_{\infty }\cdot ||\mathbf {A} ^{-1}||_{\infty }\;$, a także oznaczymy $\alpha =\epsilon K_{\infty }O\left({{9} \over {2}}n^{3}\right)\;$, a $O\left({{9} \over {2}}n^{3}\right)={{9} \over {2}}n^{3}+{{61} \over {2}}n^{2}-18n-16\;$, wtedy błąd względny wektora x jest napisany według nierówności:

${{||\delta \mathbf {x} ||_{\infty }} \over {||\mathbf {x} ||_{\infty }}}\leq {{\alpha } \over {1-\alpha }}\;$

(5.79)

Podamy teraz twierdzenie opisującą eliminację metody GCW.

Twierdzenie
Dla macierzy A, która jest nieosobliwa i jest diagonalnie dominująca kolumnowo, to w metodzie eliminacji GCW nie musimy przedstawiać wierszy.

Dowód
Załóżmy, że wyeliminowaliśmy N wierszy, a pozostałe o numerach kolumny N+1,N+2,..,n tworzą macierz diagonalnie dominująco kolumnowo. Dla macierzy $\mathbf {A} ^{(1)}\;$ poszczególne elementy macierzy L możemy napisać jako $l_{k1}=a_{ki}^{(1)}/a_{11}^{(1)}\;$ dla k=2,3,..,n, przy której z definicji macierzy diagonalnie dolinującej kolumnowo (5.47) możemy napisać warunek $\sum _{k=1}^{n}|l_{k1}|\leq 1\;$. Elementy o numerach w wierszu 2,3,..,n możemy przepisać w formie $a_{ki}^{(2)}=a_{ki}^{(1)}-l_{k1}a_{i1}^{(1)}\;$ dla numerów kolumn k=2,3,..,n. Dla macierzy A⁽²⁾ elementy macierzy spełniają warunek na diagonali $|a_{ii}^{(2)}|\geq |a_{ii}^{(1)}|-|l_{i1}|\cdot |a_{1i}^{(1)}|\;$, a elementy leżące poza diagonalą spełniają warunek $|a_{ii}^{(2)}|\leq |a_{ii}^{(1)}|+|l_{i1}|\cdot |a_{1i}^{(1)}|\;$. Z definicji macierzy, diagonalnie dominująco kolumnowo (5.47) i powyższych rozważań, możemy napisać nierówność:

$|a_{ii}^{(2)}|-\sum _{\overset {k=1}{\underset {k\neq i}{}}}^{n}|a_{ki}^{(2)}|\geq |a_{ii}^{(1)}|-|a_{1i}^{(1)}|\sum _{\overset {k=1}{\underset {k\neq i}{}}}^{n}|l_{k1}|-\sum _{\overset {k=1}{\underset {k\neq i}{}}}^{n}|a_{ki}^{(1)}|-|a_{1i}^{(1)}|\sum _{\overset {k=1}{\underset {k\neq i}{}}}^{n}|l_{k1}|=|a_{ii}^{(1)}|-\sum _{\overset {k=1}{\underset {k\neq i}{}}}^{n}|a_{ki}^{(1)}|\geq 0\;$

(5.80)

Z właśności, jeśli macierz A⁽¹⁾ jest diagonalnie dominująca kolumnowo, to macierz A⁽²⁾ jest diagonalnie dominująca kolumnowo, ogólnie rzecz biorąc macierz A^(k) dla k=1,2,..,n jest macierzą diagonalnie dominująca kolumnowo.

Rozwiązania równań liniowych metodą eliminacji Jordana (metodą eliminacji zupełnej)

Weźmy sobie układ równań, który przedstawimy w postaci macierzowej $\mathbf {A} ^{(1)}\mathbf {x} =\mathbf {b} ^{(1)}\;$, który możemy rozbić na n równań odpowiadających temu równaniu macierzowemu:

$a_{11}^{(1)}x_{1}+a_{12}^{(1)}x_{2}+...+a_{1n}^{(1)}x_{n}=b_{1}^{(1)}\;$

$a_{21}^{(1)}x_{1}+a_{22}^{(1)}x_{2}+...+a_{2n}^{(1)}x_{n}=b_{2}^{(1)}\;$
........................................................

$a_{n1}^{(1)}x_{1}+a_{n2}^{(1)}x_{2}+...+a_{nn}^{(1)}x_{n}=b_{n}^{(1)}\;$

(5.81)

Pierwsze równanie układu równań (5.81) dzielimy obustronnie przez $a_{11}^{(1)}\;$, a następnie od i-tego wiersza odejmujemy pierwszy wiersz pomnożonej przez $a_{i1}^{(1)}\;$, wtedy otrzymujemy następny układ równań, który zapisujemy w formie macierzowym $\mathbf {A} ^{(2)}\mathbf {x} =\mathbf {b} ^{(2)}\;$ rozpisując je w postaci n równań liniowych:

$x_{1}+a_{12}^{(2)}x_{2}+...+a_{1n}^{(2)}x_{n}=b_{1}^{(2)}\;$

$a_{22}^{(2)}x_{2}+...+a_{2n}^{(2)}x_{n}=b_{1}^{(2)}\;$
...........................................

$a_{n2}^{(2)}x_{2}+...+a_{nn}^{(2)}x_{n}=b_{n}^{(2)}\;$

(5.82)

Po (n-1) dokonanych eliminacjach otrzymujemy n równań algebraicznych linowych równoważnych równaniu (5.81), z których w sposób bardzo łatwy możemy policzyć niewiadome występujące w poniższym równaniu:

$x_{1}+a_{12}^{(n)}x_{3}+...+a_{1n}^{(n)}x_{n}=b_{1}^{(n)}\;$

$x_{2}+a_{23}^{(n)}x_{3}+...+a_{2n}^{(n)}x_{n}=b_{2}^{(n)}\;$
..............................................

$x_{n}=b_{n}^{(n)}\;$

(5.83)

Metoda eliminacji Jordana wymaga mnożeń $M={{1} \over {2}}n^{3}+{{1} \over {2}}n^{2}\;$ i dodawań $D={{1} \over {2}}n^{3}-{{1} \over {2}}\;$. Jak widzimy, że metoda eliminacji Jordana wymaga około półtora więcej operacji niż metoda eliminacji Gaussa. Aby w tej metodzie nie nastąpiło dzielenie przez zero należy dokonać odpowiedni wybór elementu podstawowego.

Rozkład macierzy symetrycznej A na LDL^T i LL^T

Weźmy sobie macierz kwadratową symetryczną $\mathbf {A} \;$, którą rozłożymy na iloczyn dwóch macierzy $\mathbf {L} \;$, która jest macierzą dolną trójkątną z jedynkami na diagonali i macierzy ${\overline {\mathbf {U} }}\;$, czyli: $\mathbf {L} {\overline {\mathbf {U} }}\;$, wtedy łatwo uzyskać rozkład macierzy symetrycznej w postaci $\mathbf {A} =\mathbf {L} \mathbf {D} \mathbf {U} \;$, a macierz $\mathbf {D} \;$ jest macierzą diagonalną o elementach na diagonali macierzy ${\overline {\mathbf {U} }}\;$, a $\mathbf {U} \;$ jest macierzą górną trójkątną, która jest macierzą na diagonali z jedynkami. Ponieważ $\mathbf {A} \;$ jest macierzą symetryczną, to również dobrze można napisać $\mathbf {L} \mathbf {D} \mathbf {U} =\mathbf {A} \;$, a także $\mathbf {U} ^{T}\mathbf {D} \mathbf {L} ^{T}=\mathbf {A} ^{T}=\mathbf {A} \;$, co stąd oczekujemy, że wyjdzie $\mathbf {U} =\mathbf {L} ^{T}\;$, zatem dla macierzy symetrycznych spełniony jest związek:

$\mathbf {A} =\mathbf {L} \mathbf {D} \mathbf {L} ^{T}\;$

(5.84)

Z równania (5.84) możemy powiedzieć d₁=a₁₁, a także możemy otrzymać dwa równania na l_ij i d_i dla i=2,3,..,n, i j=1,2,..,i-1 przy definicji c_ij=d_jl_ij , takiego że:

$l_{ij}=\left(a_{ij}-\sum _{k=1}^{j-1}c_{ik}l_{jk}\right)/d_{j}\;$ (5.85)

$d_{i}=a_{ii}-\sum _{k=1}^{i-1}c_{ik}l_{ik}\;$ (5.86)

W celu wyznaczenia rozkładu (5.84) należy dokonać $M={{1} \over {6}}n^{3}+n^{2}-{{7} \over {6}}n\;$ mnożeń i $D={{1} \over {6}}n^{3}-{{1} \over {6}}n\;$ dodawań. Układ równań opisanej przez równanie macierzowe $\mathbf {A} \mathbf {x} =\mathbf {b} \;$ przy definicji macierzy (5.84) możemy rozłożyć na układ trzech równań macierzowych:

$\mathbf {L} ^{T}\mathbf {z} =\mathbf {b} {\mbox{, }}\mathbf {D} \mathbf {y} =\mathbf {z} {\mbox{, }}\mathbf {L} \mathbf {x} =\mathbf {y} \;$

(5.87)

Ponieważ macierz $\mathbf {A} \;$ jest macierzą symetryczną, znając dolną część macierzy wraz z wyrazami na diagonali zużyje nam ${{1} \over {2}}n^{2}+{{1} \over {2}}n\;$ komórek pamięci komputera.

Znamy jeszcze inny rozkład macierzy $\mathbf {A} \;$, która jest zwana rozkładem Banachiewicza, która jest również zwana rozkładem Cholesky'ego:

$\mathbf {A} =\mathbf {L} \mathbf {L} ^{T}\;$

(5.88)

Zakładamy, że zachodzi tożsamość ${\overline {\mathbf {L} }}=-\mathbf {L} \;$, co z oczywistych powodów możemy napisać równoważnie dla (5.88) $\mathbf {A} ={\overline {\mathbf {L} }}{\overline {\mathbf {L} }}^{T}\;$, ale $\mathbf {L} \;$ jest macierzą trójkątną dolną niekoniecznie z jedynkami na diagonali. Wzory na l_ii i l_ij dla i=1,2,..,n i dla j=i+1,i+2,..,n piszemy:

$\sum _{k=1}^{i}l_{ij}^{2}=a_{ii}\Rightarrow l_{ii}={\sqrt {a_{ii}-\sum _{k=1}^{i-1}l_{ik}^{2}}}\;$ (5.89)

$l_{ji}=\left(a_{ji}-\sum _{k=1}^{i-1}l_{jk}l_{ik}\right)/l_{ii}\;$ (5.90)

Dla spełnionej pierwszej równości (5.89) liczby l_ij mają ograniczoną wartość dla ograniczonych elementów a_ij. Błąd zaokrągleń Banachiewicza rozkładu LL^T dla macierzy A piszemy $\mathbf {E} =\mathbf {L} \mathbf {L} ^{T}-\mathbf {A} \;$, wtedy błąd ograniczenia na tą wielkość piszemy przez:

$||\mathbf {E} ||_{2}\leq \epsilon n^{{3} \over {2}}||\mathbf {A} ||_{2}\;$

(5.91)

$||\mathbf {E} ||_{2}\leq \epsilon n^{2}||\mathbf {A} ||_{1\infty }\;$

(5.92)

Równanie macierzowe z macierzą trójdiagonalną

Weźmy sobie równanie $\mathbf {A} \mathbf {x} =\mathbf {b} \;$, w którym macierz $\mathbf {A} \;$ jest macierzą trójdiagonalną napisaną na w sposób:

$\mathbf {A} ={\begin{bmatrix}b_{1}&c_{1}&&0&\\a_{2}&b_{2}&c_{2}&&\\0&a_{3}&b_{3}&\ddots &\\&&\ddots &\ddots &c_{n-1}\\&0&&a_{n}&b_{n}\end{bmatrix}}\;$

(5.93)

Macierz $\mathbf {A} \;$ można rozłożyć na iloczyn dwóch macierzy $\mathbf {L} \;$ i $\mathbf {U} \;$, czyli na $\mathbf {L} \mathbf {U} \;$, które mają specyficzny wygląd:

$\mathbf {L} ={\begin{bmatrix}1&&0&\\l_{2}&\ddots &\\&\ddots &\ddots \\&0&l_{n}&1\end{bmatrix}}\;$

(5.94)

$\mathbf {U} ={\begin{bmatrix}u_{1}&c_{1}&&0&\\&\ddots &\ddots &\\&&\ddots &c_{n-1}\\&0&&u_{n}\\\end{bmatrix}}$

(5.95)

Wyznaczmy iloczyn macierzy $\mathbf {L} \;$ (5.94) i $\mathbf {U} \;$ (5.95), co w rezultacie po krótkich obliczeniach otrzymujemy:

$\mathbf {L} \mathbf {U} ={\begin{bmatrix}1&&0&\\l_{2}&\ddots &\\&\ddots &\ddots \\&0&l_{n}&1\end{bmatrix}}{\begin{bmatrix}u_{1}&c_{1}&&0&\\&\ddots &\ddots &\\&&\ddots &c_{n-1}\\&0&&u_{n}\\\end{bmatrix}}={\begin{bmatrix}u_{1}&c_{1}&\\l_{2}u_{1}&l_{2}c_{1}+u_{2}&c_{2}\\&\ddots &\ddots &c_{n-1}\\&&l_{n}u_{n-1}&l_{n}c_{n-1}+u_{n}\end{bmatrix}}$

(5.96)

Patrząc na obliczenia (5.96) i porównując ten wynik z (5.93), to możemy podać ogólny wynik na l_i i u_i, który podajemy w postaci przepisu:

$u_{1}=b_{1}\;$

(5.97)

$l_{i}={{a_{i}} \over {u_{i-1}}}\;$

(5.98)

$u_{i}=b_{i}-l_{i}c_{i-1}{\mbox{ dla }}i=2,3,..,n\;$

(5.99)

Mamy sobie równanie macierzowe $\mathbf {A} \mathbf {x} =\mathbf {d} \;$, który po rozkładzie macierzy $\mathbf {A} \;$ na iloczyn dwóch czynników macierzowych (5.94) i (5.95) zapisując jako równanie $\mathbf {L} \mathbf {U} \mathbf {x} =\mathbf {d} \;$, który rozbijemy na dwa równania, tzn. na $\mathbf {L} \mathbf {y} =\mathbf {d} \;$ i $\mathbf {U} \mathbf {x} =\mathbf {y} \;$, z których pierwsze równanie rozpiszemy w postaci równania y_i jako:

$y_{1}=d_{1}\;$

(5.100)

$y_{i}=d_{i}-l_{i}y_{i-1}{\mbox{ dla }}i=2,3,..,n\;$

(5.101)

a drugie równanie macierzowe rozpisujemy na element x_i w zależności od zmiennej y_i i x_i+1 i u_i:

$x_{n}={{y_{n}} \over {u_{n}}}\;$

(5.102)

$x_{i}=(y_{i}-c_{i}x_{i+1})/u_{i}\;$

(5.103)

Gdy $\mathbf {A} \;$ jest macierzą diagonalnie dominująca kolumnowo, tzn. gdy zachodzą warunki: |b₁|≥|a₂|, |b_i|≥|c_i-1|+|a_i+1| dla i=2, 3,..,n-1, |b_n|≥|c_n-1|, to ten nasz rozkład jest oparty z częściowych wyborem elementu podstawowego, a więc jest metodą niezawodną, tzn. przy niej nie występuje dzielenie przez zero. Natomiast, gdy macierz $\mathbf {A} \;$ jest macierzą diagonalnie dominująca, tzn. gdy zachodzą nierówności: |b₁|≥|c₁|, |b_i|≥|a_i|+|c_i| dla i=2,3,..,n-1, |b_n|≥|a_n|, to również metoda eliminacji Gaussa w tym przypadku jest niezawodna.

Weźmy sobie inny rozkład macierzy $\mathbf {A} \;$ na iloczyn dwóch czynników, które te macierze są zdefiniowane w formie:

$\mathbf {L} ={\begin{bmatrix}l_{1}&&0&\\a_{2}&l_{2}\\&\ddots &\ddots \\&0&a_{n}&l_{n}\end{bmatrix}}\;$

(5.104)

$\mathbf {U} ={\begin{bmatrix}1&u_{1}&&0&\\&1&u_{2}\\&&\ddots &\ddots \\&&&&u_{n-1}\\&0&&&1\end{bmatrix}}\;$

(5.105)

Policzmy teraz iloczyn macierzy $\mathbf {L} \;$ (5.104) i $\mathbf {U} \;$ (5.105), który napiszemy według przepisu:

$\mathbf {L} \mathbf {U} ={\begin{bmatrix}l_{1}&&0&\\a_{2}&l_{2}\\&\ddots &\ddots \\&0&a_{n}&l_{n}\end{bmatrix}}{\begin{bmatrix}1&u_{1}&&0&\\&1&u_{2}\\&&\ddots &\ddots \\&&&&u_{n-1}\\&0&&&1\end{bmatrix}}=\;$
$={\begin{bmatrix}l_{1}\\a_{2}&a_{2}u_{1}+l_{2}&l_{2}u_{2}\\\ddots &\ddots &\ddots \\&&a_{n-1}&a_{n-1}u_{n-2}+l_{n-1}&l_{n-1}u_{n-1}\\&0&&a_{n}&a_{n}u_{n-1}+l_{n}\\\end{bmatrix}}\;$

(5.106)

Patrząc na obliczenia (5.106) i porównując ten wynik z (5.93), to możemy podać ogólny wynik na l_i i u_i, który podajemy w postaci przepisu:

$l_{1}=b_{1}\;$

(5.107)

$u_{1}={{c_{i}} \over {l_{i}}}\;$

(5.108)

$l_{i}=b_{i}-a_{i}u_{i-1}\;$

(5.109)

Mamy sobie równanie macierzowe $\mathbf {A} \mathbf {x} =\mathbf {d} \;$, który po rozkładzie macierzy $\mathbf {A} \;$ na iloczyn dwóch czynników macierzowych (5.104) i (5.105) zapisując jako równanie $\mathbf {L} \mathbf {U} \mathbf {x} =\mathbf {d} \;$, który rozbijemy na dwa równania, tzn. na $\mathbf {L} \mathbf {y} =\mathbf {d} \;$ i $\mathbf {U} \mathbf {x} =\mathbf {y} \;$, z których pierwsze równanie rozpiszemy w postaci równania y_i dla i=1,2,..,n-1 jako:

$y_{1}={{d_{1}} \over {l_{1}}}\;$

(5.110)

$y_{i+1}=(d_{i+1}-a_{i+1})/l_{i+1}\;$

(5.111)

a drugie równanie macierzowe rozpisujemy na element x_i w zależności od zmiennej y_i i x_i+1 i u_i dla i=n-1,n-2,..,1:

$x_{n}=y_{n}\;$

(5.112)

$x_{i}=y_{i}-u_{i}x_{i+1}\;$

(5.113)

Będziemy szukali błędów zaokrągleń przy liczeniu macierzy $\mathbf {A} \;$, którą liczymy ze wzoru $\mathbf {A} =\mathbf {L} \mathbf {U} \;$. Błąd zaokrągleń będziemy liczyli w naszym przypadku ze wzoru $\mathbf {E} =\mathbf {L} \mathbf {U} -\mathbf {A} \;$, wtedy dla p=1,∞ możemy napisać ||E||_p≤2ε||T||_p. Jeśli natomiast (L+δL)y=d i (U+δU)x=y, wtedy powiemy: ||δL||_p≤3ε, ||L||_p≤2ε, a także zachodzą ||δU||_p≤3ε||T||_p i ||U||_p≤2ε||T||_p. Jeśli oznaczymy (A+δA)x=d, i ||δT||_p≤14ε||T||_p przy definicji ||A^-1||_p||A||_p=K_p dla p=1,∞, a także wprowadzimy oznaczenie na parametr α=14εK_p, wtedy wzór (5.36) przyjmuje taki sam wygląd jak wzór (5.79).

Równanie macierzowe z macierzą podobną do trójdiagonalnej

W matematyce spotykamy się też w równaniach macierzowych $\mathbf {A} \mathbf {x} =\mathbf {b} \;$ z macierzami $\mathbf {A} \;$ podobnej do macierzy trójdiagonalnej (5.94), której przepis jest:

$\mathbf {A} ={\begin{bmatrix}b_{1}&c_{1}&\cdots &\cdots &\cdots &p_{1}\\a_{2}&b_{2}&c_{2}&&0&\vdots \\\vdots &a_{3}&b_{3}&\ddots &&\vdots \\\vdots &&\ddots &\ddots &\ddots &\vdots \\\vdots &0&&\ddots &\ddots &c_{n-1}\\q_{1}&\cdots &\cdots &\cdots &a_{n}&b_{n}\\\end{bmatrix}}\;$

(5.114)

W porównaniu z macierzą trójdiaagonalną macierz powyższa róźni się od poprzedniej tym, że powyżej pojawiają się niezerowe elementy p₁ i q₁. Rozłóżmy macierz (5.114) na iloczyn $\mathbf {L} \mathbf {U} \;$ stosując metodę Gaussa lub Doolittle'a. Macierz $\mathbf {L} \;$ jest macierzą trójkatną dolną z jedynkami na diagonali, a macierz $\mathbf {U} \;$ jest macierzą trójkątną górną niekoniecznie z jedynkami na diagonali, które przedstawiamy w formie:

$\mathbf {L} ={\begin{bmatrix}1&\\l_{2}&\ddots &&0&&\\&\ddots &\ddots &\\&0&l_{n-1}&\ddots &\\q_{2}&\cdots &q_{n-1}&q_{n}&1\end{bmatrix}}\;$

(5.115)

$\mathbf {U} ={\begin{bmatrix}u_{1}&c_{1}&\cdots &\cdots &p_{1}\\&\ddots &\ddots &0&\vdots \\&&\ddots &c_{n-1}&p_{n-2}\\&0&&\ddots &p_{n-1}\\&&&&u_{n}\end{bmatrix}}$

(5.116)

Policzmy teraz iloczyn macierzy $\mathbf {L} \;$ (5.115) i $\mathbf {U} \;$ (5.116) i przyrównamy je do macierzy podobnej do trójdiagonalnej $\mathbf {A} \;$ (5.114):

$\mathbf {L} \mathbf {U} ={\begin{bmatrix}1&\\l_{2}&\ddots &&0&&\\&\ddots &\ddots &\\&0&l_{n-1}&\ddots &\\q_{2}&\cdots &q_{n-1}&q_{n}&1\end{bmatrix}}{\begin{bmatrix}u_{1}&c_{1}&\cdots &\cdots &p_{1}\\&\ddots &\ddots &0&\vdots \\&&\ddots &c_{n-1}&p_{n-2}\\&0&&\ddots &p_{n-1}\\&&&&u_{n}\end{bmatrix}}=\;$
$={\begin{bmatrix}u_{1}&c_{1}&&&&p_{1}\\u_{1}l_{2}&l_{2}c_{1}+u_{2}&c_{2}&&&\vdots \\\vdots &u_{2}l_{3}&l_{3}c_{2}+u_{3}&\ddots &&\vdots \\\vdots &&\ddots &\ddots &\ddots &\vdots \\\vdots &0&&\ddots &\ddots &p_{n-1}+l_{n-1}p_{n-2}\\u_{1}q_{2}&q_{2}c_{1}+q_{3}u_{2}&q_{3}c_{2}+q_{4}u_{3}&&&q_{2}p_{1}+\dots +p_{n-1}q_{n}+...+u_{n}\\\end{bmatrix}}$

(5.117)

Porównując macierz końcową uzyskaną w wyniku obliczeń (5.117) i macierz (5.114), wtedy możemy napisać tożsamości na parametry p_i, u_i i q_i:

Twierdzenie Wzory na elementy macierzy L i U

$u_{i}=b_{1}\;$

$l_{2}=a_{2}/u_{1}\;$

$q_{2}=-q_{1}/u_{1}\;$

---

$u_{2}=b_{2}-l_{2}c_{1}\;$

$l_{3}=a_{3}/u_{2}\;$

$q_{3}=-q_{2}c_{1}/u_{2}\;$

---

$u_{n-2}=b_{n-2}-l_{n-2}c_{n-3}\;$

$l_{n-1}=a_{n-1}/u_{n-2}\;$

$q_{n-1}=-q_{n-2}c_{n-3}/u_{n-2}\;$

---

$u_{n-1}=b_{n-1}-l_{n-1}c_{n-2}\;$

$q_{n}=(a_{n}-q_{n-1}c_{n-2})/u_{n-1}\;$

---

$p_{2}=-l_{2}p_{1}\;$

$p_{3}=-l_{3}p_{2}\;$

---

$p_{n-2}=-l_{n-2}p_{n-3}\;$

$p_{n-1}=c_{n-1}-l_{n-1}p_{n-2}\;$

$u_{n}=b_{n}-p_{1}q_{2}-...-p_{n-1}q_{n}\;$

Mamy sobie równanie macierzowe $\mathbf {A} \mathbf {x} =\mathbf {d} \;$, który po rozkładzie macierzy $\mathbf {A} \;$ na iloczyn dwóch czynników macierzowych (5.115) i (5.116) zapisujemy je jako równanie $\mathbf {L} \mathbf {U} \mathbf {x} =\mathbf {d} \;$, który rozbijemy na dwa równania, tzn. na $\mathbf {L} \mathbf {y} =\mathbf {d} \;$ i $\mathbf {U} \mathbf {x} =\mathbf {y} \;$, z których pierwsze równanie rozpiszemy w postaci równania y_i, a drugie x_i dla i=1,2,..,n:

$y_{1}=d_{1}\;$
$y_{2}=d_{2}-l_{2}y_{1}\;$
$\cdots \cdots \cdots \;$
$y_{n-1}=d_{n-1}-l_{n-1}y_{n-2}\;$
$y_{n}=d_{n}-q_{2}y_{1}-...-q_{n}y_{n-1}\;$

(5.118)

$x_{n}=y_{n}/u_{n}\;$
$x_{n-1}=(y_{n-1}-p_{n-1}x_{n})/u_{n-1}\;$
$x_{n-2}=(y_{n-2}-c_{n-2}x_{n-1}-p_{n-2}x_{n})/u_{n-2}\;$
$\cdots \cdots \cdots \;$
$x_{1}=(y_{1}-c_{1}x_{2}-p_{1}x_{n})/u_{1}\;$

(5.119)

Podczas zaokrąglania macierzy przy liczeniu iloczynu LU powstaje błąd zaokrąglenia E=LU-A, wtedy norma błędu jest opisywana przez ||E||₁≤ 2nε||A||₁. Jeśli wyniki L i U są w pewien w sposób zaokrąglone, wtedy możemy napisać rozwiązanie równań macierzowych (L+δL)/y=d i (U+δU)/x=y, wtedy ograniczenia na normy i normy błędów dla macierzy L i U, przy zaokrągleniu, piszemy ||L||₁≤2, ||δL||₁≤ε||L||₁, a także ||U||₁≤2||A||₁, ||δU||₁≤5ε||U||₁. Norma błędu zaokrąglenia macierzy A przedstawiamy w zależności od normy tejże macierzy w sposób: ||δA||₁≤(6n+20+20nε)||A||₁. Patrząc na wynik na normą błędu zaokrąglenia ostatnio napisany w tym rozdziale dowiadujemy się że nastąpiło pogorszenie osaczowania obliczeń dla naszej macierzy A w tym rozdziale niż dla macierzy całej A według jego normy błędu (5.78).

Wyznaczanie wartości wyznacznika oraz macierzy odwrotnej

Podczas wyznaczania wyznacznika macierzy należy wykonać n! obliczeń, co jest za dużo dla obecnych maszyn cyfrowych i dlatego stosuje się rozkład macierzy A na dwa macierze L i U do postaci A=LU według metody eliminacji Gaussa, wtedy po tym rozkładzie zastosujmy twierdzenie o wyznaczniku iloczynu dwóch macierzy:

$\operatorname {det} (\mathbf {A} +\mathbf {E} )=\operatorname {det} \mathbf {L} \operatorname {det} \mathbf {U} =\operatorname {det} \mathbf {U} \;$

(5.120)

Powyżej zastosowaliśmy to , że detL=1, oraz że macierz LU została policzona z błędem na maszynie cyfrowej z błędem E. Przybliżeniu wyznacznik macierzy A w przybliżeniu jest równy iloczynowi n elementów na przekątnej macierzy U.

Oznaczmy przez X odwrotność macierzy A, wtedy błąd liczenia iloczynu XA powinien być w przybliżeniu równy macierzy jednostkowej z błędem określonej E=I-XA, wtedy błąd bezwzględny przy wyznaczaniu X według normy ||⋅||_∞ jest równy:

${{||\mathbf {E} ||_{\infty }} \over {||\mathbf {X} ||_{\infty }}}\leq \epsilon g\left(2,005n^{2}+n^{3}+{{1} \over {4}}n^{4}\epsilon \right)\;$

(5.121)

przez g oznaczymy największy element układu macierzy A⁽¹⁾,A⁽²⁾,...,A⁽ⁿ⁾=U, które są otrzymane metodą eliminacji jakie dotychczas poznaliśmy. Wartość g możemy napisać według oszacowania g≤8||A||_1∞. Jeśli mamy oszacowanie odwrotności macierzy A, wtedy błąd zmiennej x przepisujemy według (5.79). Oszacowanie błędu E według normy ||⋅||₂ piszemy $||\mathbf {E} ||_{2}\leq \epsilon n^{{3} \over {2}}||\mathbf {A} ||_{2}\;$. Jeśli dokonamy rozkładu macierzy A według rozkładu Banachiewicza, to oznaczmy wartości własne macierzy A przez symbole λ₁,...,λ_n. a wartości własne macierzy LL^T napiszemy przez γ₁,...,γ_n. Napiszmy przez ||A||₂=λ₁, ||LL^T||₂=γ₁ jako największe wartości własne wspomnianych wyżej macierzy, wtedy na podstawie wyżej wniosków:

$-\epsilon n^{{3} \over {2}}\lambda _{1}\leq \lambda _{i}-\gamma _{i}\leq \epsilon n^{{3} \over {2}}\lambda _{1}\;$

(5.122)

Z nierówności (5.122) możemy napisać następną dalszą nierówność, który piszemy względem wyznacznika macierzy LL^T, wtedy:

$\prod _{i=1}^{n}\left(\lambda _{i}-\epsilon n^{{3} \over {2}}\lambda _{1}\right)\leq \operatorname {det} (\mathbf {L} \mathbf {L} ^{T})\leq \prod _{i=1}^{n}\left(\lambda _{i}-\epsilon n^{{3} \over {2}}\lambda _{1}\right)\Rightarrow \;$
$\Rightarrow \prod _{i=1}^{n}\left(1-\epsilon n^{{3} \over {2}}{{\lambda _{1}} \over {\lambda _{i}}}\right)\leq {{\operatorname {det} (\mathbf {L} \mathbf {L} ^{T})} \over {\operatorname {det} \mathbf {A} }}\leq \prod _{i=1}^{n}\left(1-\epsilon n^{{3} \over {2}}{{\lambda _{1}} \over {\lambda _{i}}}\right)$

(5.123)

Na podstawie obliczeń (5.123) błąd bezwzględny wyznaczania wyznacznika det(LL^T) jest opisywany poprzez wyrażenie εn^3/2λ_max/λ_min=εn^3/2K₂. Podobnie postępujemy przy wyznaczaniu błędu bezwzględnego wyznacznika A^-1 według normy ||⋅||₂.

${{||\mathbf {A} ^{-1}-(\mathbf {L} ^{T})^{-1}(\mathbf {L} ^{-1})||_{2}} \over {||\mathbf {|} |\mathbf {A} ^{-1}||_{2}}}\leq \epsilon n^{{3} \over {2}}K_{2}\;$

(5.124)

Przy wyznaczaniu odwrotności macierzy A, którą możemy rozłożyć na L i U, wtedy macierz odwrotną macierzy możemy policzyć ze wzoru LUx⁽ⁱ⁾=e⁽ⁱ⁾, gdzie e⁽ⁱ⁾ jest macierzą kanoniczną, gdzie na i-tym miejscu jest jedynka, a na pozostałych elementach są zera. Element x⁽ⁱ⁾ jest i-tą kolumną macierzy odwrotnej do LU. Przy liczeniu macierzy odwrotnej (5.94) należy skorzystać ze wzoru Lx⁽ⁱ⁾=e⁽ⁱ⁾ dla i=1,2,..,n, wtedy błąd względny przy liczeniu błędu macierzy odwrotnej do L przy normie ||⋅||₁ opisujemy przez:

${{||\operatorname {fl} (\mathbf {L} ^{-1})-\mathbf {L} ^{-1}||_{1}} \over {||\mathbf {L} ^{-1}||_{1}}}\leq \epsilon (n+2)K_{1}\;$

(5.125)

Odwrotność macierzy diagonalnej do (5.94) jest macierzą trójkątną dolną z jedynkami na diagonali zapisaną:

$\mathbf {L} ^{-1}={\begin{bmatrix}1&&&&0&\\-l_{2}&1\\l_{2}l_{3}&-l_{3}&\ddots \\\vdots &\ddots &\ddots \\(-1)^{n-1}l_{2}...l_{n}&(-1)^{n-2}l_{3}...l_{n}&\cdots &-l_{n-1}l_{n}&-l_{n}&1\\\end{bmatrix}}\;$

(5.126)

Aby wyznaczyć wektor y=L^-1x należy wykonać n(n-1)/2 mnożeń, a przy liczeniu wektora Ly=x należy wykonać n-1 mnożeń.

Poprawianie rozwiązań układów równań liniowych i wektor reszt

Podczas rozwiązywania układów równań liniowych Ax=b jest na ogół obarczona ściśle określonym pewnym błędem. W celu napisania jaki błąd popełniliśmy należy obliczyć resztę:

$\mathbf {r} =\mathbf {b} -\mathbf {A} \mathbf {x} \;$

(5.127)

przy którym sprawdzimy, czy ona jest równa zeru. Przy metodach dokładnych podczas dokonanych przybliżeń, w tychże metodach ta reszta jest rzeczywiście jest nierówna zero. Dla rozwiązania x równania macierzowego liniowego Ax=b potrafimy dokładnie wyznaczyć poprawkę do uzyskanego rozwiązania w wyniku zaokrągleń, wyniku czego dokładne rozwiązanie zapisujemy $\mathbf {x} +\delta \mathbf {x} ={\hat {\mathbf {x} }}\;$, co dokładną poprawką δx jest rozwiązaniem równania Aδx=r. Dla poszukiwanego rozwiązania metodą GCW uzyskując przedtem macierze L i U wykonując przy tym n² mnożeń i n²-n dodawań, ale wtedy uzyskamy przybliżoną wartość błędu δx+δ(δx), zatem dla poprawionego rozwiązania naszego równania macierzowego liniowego ${\overline {\mathbf {x} }}=\mathbf {x} +\delta \mathbf {x} +\delta (\delta \mathbf {x} )\;$ wektor reszt będzie miał normę z oczywistych powodów mniejszą wartość. Sformujemy twierdzenie, który coś mówi o normie ||⋅||_∞; reszty z dzielenia względem wielkości ε:

Lemat
Jeśli dla wektora reszty obliczonej dokładnie wyznaczonej według równości (5.127) poprawka do rozwiązania przybliżonego δx została wyznaczona metodą GCW, dla której zachodzi nierówność:

${{1} \over {2}}||\mathbf {r} ||_{\infty }\geq \epsilon ||\mathbf {A} ||_{\infty }\left({{9} \over {2}}n^{3}+{{61} \over {2}}n^{2}-18n-16\right)||\delta \mathbf {x} ||_{\infty }+||\mathbf {A} \mathbf {x} ||_{\infty }\epsilon \;$

(5.128)

wtedy jest spełniony warunek względem poprawionego rozwiązania ${\overline {\mathbf {x} }}\;$ i rozwiązania przybliżonego x:

$||\mathbf {b} -\mathbf {A} {\overline {x}}||_{\infty }\leq {{1} \over {2}}||\mathbf {b} -\mathbf {A} \mathbf {x} ||\;$

(5.129)

Macierzowe algebraiczne liniowe równania iteracyjne

Określmy sobie ciąg wektorów x⁽⁰⁾,x⁽¹⁾,...,x⁽ⁱ⁾, które są określone przez równość iteracyjną elementu o numerze i+1 w zależności od elementu i-tego dla i=0,1,... według przepisu:

$\mathbf {x} ^{(i+1)}=\mathbf {M} \mathbf {x} ^{(i)}+\mathbf {w} \;$

(5.130)

Twierdzenie
Jeśli mamy równanie iteracyjne (5.130), to on jest zbieżny przy ρ(M)<1.

Dowód
Według równania iteracyjnego (5.130) napiszemy wzór na i+1-ty element zmiennej x⁽ⁱ⁾, wtedy napiszemy:

$\mathbf {x} ^{(i)}=\mathbf {M} \mathbf {x} ^{(i)}+\mathbf {w} =\mathbf {M} (\mathbf {M} +\mathbf {w} )+\mathbf {w} =\mathbf {M} ^{(i+1)}\mathbf {x} ^{(0)}+\left(\mathbf {M} ^{(i)}\mathbf {w} +\mathbf {M} ^{(i-1)}\mathbf {w} +...+\mathbf {w} \right)\;$

(5.131)

Ponieważ ρ(M)<1 i $\mathbf {M} ^{(i+1)}\mathbf {x} ^{(0)}{\xrightarrow[{i\rightarrow \infty }]{}}0\;$, wtedy szereg $\mathbf {M} ^{(i)}\mathbf {w} +\mathbf {M} ^{(i-1)}\mathbf {w} +...+\mathbf {w} \;$ jest zbieżny do pewnego szeregu, co stąd dla tego określonego M szereg generowany przez równanie iteracyjne (5.130) jest szeregiem zbieżnym.

Według równania (5.131) przy zachodzącym ρ(M)<1 szereg jest szeregiem zbieżnym do x^(∞) przy zachodzącym równaniu x^(∞)=Mx^(∞)+w. Przy zachodzącym rozwiązywaniu równania AX=b należy tak dobrać takie M by był spełniony warunek w powyższym twierdzeniu. Ponieważ zachodzi x=Mx+w, wtedy wyznaczając x z równania liniowego względem wektora zmiennych i podstawiając do granicznego równania uzyskanego z (5.130), otrzymujemy:

$\mathbf {x} =\mathbf {M} \mathbf {x} +\mathbf {w} \Rightarrow (\mathbf {I} -\mathbf {M} )\mathbf {x} =\mathbf {w} \Rightarrow \mathbf {w} =(\mathbf {I} -\mathbf {M} )\mathbf {A} ^{-1}\mathbf {b} \;$

(5.132)

Przyjmijmy tożsamość macierzową M=I-NA, wtedy tożsamość macierzową (5.132) możemy podstawić do równania iteracyjnego (5.130) otrzymując następną, ale równoważną tożsamość iteracyjną:

$\mathbf {x} ^{(i+1)}=(\mathbf {I} -\mathbf {N} \mathbf {A} )\mathbf {x} ^{(i)}+\mathbf {N} \mathbf {b} \;$

(5.133)

W każdej iteracji Mx⁽ⁱ⁾+w wyliczamy x⁽ⁱ⁺¹⁾=Mx⁽ⁱ⁾+w+δ⁽ⁱ⁾, gdzie δ⁽ⁱ⁾ jest błędem zaokrągleń, która jest wielkością bardzo małą, ale przy większej liczbie iteracji zaczyna odgrywać ogromną rolę:

$\mathbf {x} ^{(i+1)}=\mathbf {M} ^{(i+1)}\mathbf {x} ^{(0)}+\mathbf {M} ^{(i)}\mathbf {w} +...+\mathbf {w} +\delta ^{(i)}+\mathbf {M} \delta ^{(i-1)}+...+\mathbf {M} ^{i}\delta ^{(0)}\;$

(5.134)

Błąd zaokrągleń δ⁽ⁱ⁾+Mδ^(i-1)+...+Mⁱδ⁽⁰⁾ może spowodować generowanie ciągu cyklicznego, który krótko opiszemy jako x⁽ⁱ⁺¹⁾=x⁽⁰⁾, który nie jest zbieżny do żadnego rozwiązania, przed którą sytuacją jest trudno ustrzec.

Rozwiązanie algebraicznych układów równań metodą Jacobiego

Weźmy sobie macierz $\mathbf {A} \;$ układów równań algebraicznych zapisanych w postaci macierzowej $\mathbf {A} \mathbf {x} =\mathbf {b} \;$, w której wspomnianą macierz możemy rozłożyć na trzy części $\mathbf {A} =\mathbf {L} +\mathbf {D} +\mathbf {U} \;$, w której macierz $\mathbf {L} \;$ jest macierzą poddiagonalną, $\mathbf {D} \;$ jest macierzą diagonalną, a $\mathbf {U} \;$ macierzą naddiagonalną, i przyjmując jednocześnie w (5.133), że $\mathbf {N} =\mathbf {D} ^{-1}\;$, wtedy wspomniane wyrażenie dla i=0,1,2,... możemy zapisać:

$\mathbf {x} ^{(i+1)}=\left(\mathbf {I} -\mathbf {D} ^{-1}\left(\mathbf {L} +\mathbf {D} +\mathbf {U} \right)\right)\mathbf {x} ^{(i)}+\mathbf {D} ^{-1}\mathbf {b} \Rightarrow \mathbf {D} \mathbf {x} ^{(i+1)}=-\left(\mathbf {L} +\mathbf {U} \right)\mathbf {x} ^{(i)}+\mathbf {b} \;$

(5.135)

Jeśli chcemy zastosować wzór (5.135), to równania algebraiczne równania macierzowego $\mathbf {A} \mathbf {x} =\mathbf {b} \;$ powinny mieć na diagonali elementy tylko niezerowe, a jeśli są jakieś elementy niezerowe, to wybieramy jakąś kolumnę w której znajduje się największa liczba zer, i tak przedstawiamy wiersze w tym układzie równań by element o maksymalnym module był niezerowy, by on potem po przedstawieniu znajdował się na diagonali, by potem ten wiersz pominąć, tzn. nie rozważać go w dalszych przestawianiach. Czynność tą powtarzamy dla pozostałych kolumn, w ten sposób po tych operacjach otrzymujemy na diagonali tylko niezerowe elementy. Ta metoda jest zawsze niezawodna, ale pracochłonna. Jeśli w każdej kolumnie powyżej rozważanej macierzy znajduje się taka sama liczba zerowych elementów, co mamy doczynienia z macierzami rzadkimi, to wybór elementu podstawowego postępujemy zgodnie z metodą GCW. W powyższej metodzie oczywiście mamy:

$\mathbf {M} _{J}=\mathbf {M} =\mathbf {I} -\mathbf {N} \mathbf {A} =\mathbf {I} -\mathbf {D} ^{-1}(\mathbf {L} +\mathbf {D} +\mathbf {U} )=-\mathbf {D} ^{-1}\left(\mathbf {L} +\mathbf {U} \right)\;$

(5.136)

Należy zauważyć, że zastosowanie powyższej metody nie gwarantuje zbieżności metody, tzn. ρ(D^-1(L+U))<1, ale gdy macierz A jest macierzą silnie diagonalnie dominującą lub silnie diagonalnie dominującą kolumnowo, to metoda Jacobiego jest na pewno zbieżna.

Rozwiązanie algebraicznych układów równań metodą Gaussa-Seidla

Podobnie jak w metodzie, dla takiego samego typu macierzy mamy to samo $\mathbf {A} =\mathbf {L} +\mathbf {D} +\mathbf {U} \;$, wtedy dla równania (5.133) możemy przyjąć $\mathbf {N} =\left(\mathbf {D} +\mathbf {L} \right)^{-1}\;$, wtedy wspomniane równanie możemy zapisać:

$\mathbf {x} ^{(i+1)}=\left(\mathbf {I} -\left(\mathbf {D} +\mathbf {L} \right)^{-1}\left(\mathbf {L} +\mathbf {D} +\mathbf {U} \right)\right)\mathbf {x} ^{(i)}+(\mathbf {D} +\mathbf {L} )^{-1}\mathbf {b} \Rightarrow \;$
$\Rightarrow (\mathbf {D} +\mathbf {L} )\mathbf {x} ^{(i+1)}=\left(\mathbf {D} +\mathbf {L} -\left(\mathbf {L} +\mathbf {D} +\mathbf {U} \right)\right)\mathbf {x} ^{(i)}+\mathbf {b} \Rightarrow (\mathbf {D} +\mathbf {L} )\mathbf {x} ^{(i+1)}=-\mathbf {U} \mathbf {x} ^{(i)}+\mathbf {b} \Rightarrow \;$
$\Rightarrow \mathbf {D} \mathbf {x} ^{(i+1)}=-\mathbf {L} \mathbf {x} ^{(i+1)}-\mathbf {U} \mathbf {x} ^{(i)}+\mathbf {b}$

(5.137)

Przy końcowej równości (5.136) powstaje pytanie w jaki sposób można obliczyć prawą stronę wspomnianej równości nie znając $\mathbf {x} ^{(i+1)}\;$, otóż problem jest taki, gdy chcemy obliczyć $x_{1}^{(i+1)}\;$ już nie trzeba znać innych elementów elementów tego typu, już dla $x_{j}^{(i+1)}\;$ będziemy mogli wykorzystać już policzone elementy dla numerów składowych $\mathbf {x} ^{(i+1)}\;$ mniejszych niż "j". By móc zastosować metodę Gaussa-Seidla należy pamiętać, by na diagonali znajdowały się niezerowe elementy. W powyższej metodzie oczywiście mamy:

$\mathbf {M} _{GS}=\mathbf {M} =\mathbf {I} -\mathbf {N} \mathbf {A} =\mathbf {I} -\left(\mathbf {D} +\mathbf {L} \right)^{-1}\left(\mathbf {L} +\mathbf {D} +\mathbf {U} \right)=-\left(\mathbf {D} +\mathbf {L} \right)^{-1}\mathbf {U} \;$

(5.138)

Gdy 0<ρ(M_j)<1 metoda Gaussa-Seidla jest bardziej zbieżna niż metoda Jacobiego, bo zachodzi ρ(M_GS)<ρ(M_J). Macierz $\mathbf {A} \;$ jest macierzą symetryczną, to metoda Gaussa-Seidla jest zbieżna, gdy wspomniana macierz jest dodatnio określona. Gdy macierz $\mathbf {A} \;$ jest macierzą silnie dominująca diagonalnie lub silnie dominująca diagonalnie kolumnowo metoda Gaussa-Seidla jest metodą na pewno zbieżna, ale silniej niż metoda Jacobiego.

Błędy iteracyjne w algebraicznych równaniach macierzowych

Jeśli do wzoru (5.134), który przedstawia iteracje zmiennej o wskaźniku i-tym, który jest napisany w zależności od błędów zaokrągleń podstawimy wzór (5.131), który przedstawia wzór w których obliczenia są dokonane bardzo dokładne, co w rezultacie możemy powiedzieć dla j=0,1,2,...,i:

${\tilde {\mathbf {x} }}^{(i+1)}-{\tilde {\mathbf {x} }}^{(i+1)}=\delta ^{(i)}+\mathbf {M} \delta ^{(i-1)}+...+\mathbf {M} ^{(i)}\delta ^{(0)}\;$

(5.139)

W powyższy wzór jest napisany dla ciągu wektorów x⁽¹⁾, x⁽²⁾,..., który jest zbieżny do dokładnego rozwiązania ${\hat {x}}\;$, który jest rozwiązaniem równania macierzowego $\mathbf {A} \mathbf {x} =\mathbf {b} \;$, dalej możemy przyjąć:

${{1} \over {2}}||{\tilde {\mathbf {x} }}||<||\mathbf {x} ^{(i)}||\leq ||2{\tilde {\mathbf {x} }}||{\mbox{ dla }}{\tilde {\mathbf {x} }}\neq 0$

(5.140)

wtedy możemy oszacowanie wynikające z (5.139) przy założeniu, że zachodzi ||δ⁽ⁱ⁾||<χ, napisać według:

$||{\tilde {\mathbf {x} }}^{(i+1)}-\mathbf {x} ^{(i+1)}||\leq \delta ^{(i)}+||\mathbf {M} ||\delta ^{(i-1)}+||\mathbf {M} ^{(2)}||\delta ^{(i-2)}+...+||\mathbf {M} ^{(i)}||\delta ^{(0)}\leq \;$
$\leq \left(1+||\mathbf {M} ||+||\mathbf {M} ^{(2)}||+...+||\mathbf {M} ^{(i)}||\right)\chi \Rightarrow \;$
$\Rightarrow {{||{\tilde {\mathbf {x} }}^{(i+1)}-\mathbf {x} ^{(i+1)}||} \over {||{\tilde {\mathbf {x} }}^{(i+1)}||}}\leq \left(1+||\mathbf {M} ||+||\mathbf {M} ^{(2)}||+...+||\mathbf {M} ^{(i)}||\right){{2\chi } \over {||{\tilde {\mathbf {x} }}||}}\Rightarrow \;$
$\Rightarrow {{||{\tilde {\mathbf {x} }}^{(i+1)}-\mathbf {x} ^{(i+1)}||} \over {||{\tilde {\mathbf {x} }}^{(i+1)}||}}\leq {{1} \over {1-||\mathbf {M} ||}}{{2\chi } \over {||{\tilde {\mathbf {x} }}||}}$

(5.141)

Rozwiązanie algebraicznych układów równań metodą Czebyszewa

Mamy sobie metodę iteracyjną (5.133) za pomocą której generujemy ciąg wektorów x⁽¹⁾,x⁽²⁾,..., który jest zbieżny do dokładnego rozwiązania układu równań zapisanych w sposób macierzowy Ax=b przy tak dobranym N=N_s, w taki sposób by moduł z M zapisanej jako ρ(I-N_sA)<1. Można tak założyć, że macierz A jest macierzą symetryczną i dodatnio określoną, i można to uzyskać biorąc wielomian macierzowy N_s, który powstaje po zastąpieniu w miejsce t macierzy A, który jest stopnia s-1 wówczas zbudujmy wielomian M_s, który formujemy według zasady p_s=1-w(t)t gdzie w miejsce t podstawiamy macierz A. Sformułujmy taki wielomian w(t) taki by był spełniony warunek $\max _{t\in \langle \alpha ,\beta \rangle }\left|1-w(t)t\right|\;$, by było najmniejsze jak tylko możliwe. Wielomian p_s(t), który przedstawiliśmy powyżej możemy przepisać jako:

$p_{s}(t)={{T_{s}\left({{2t-(\beta -\alpha )} \over {\beta -\alpha }}\right)} \over {T_{s}\left(-{{\beta +\alpha } \over {\beta -\alpha }}\right)}}\;$

(5.142)

Wielomiany Czybyszewa są zdefiniowane wzorem (1.27), wtedy dla s=1 wyrażenie (5.142) możemy na podstawie definicji wyżej wspomnianych wielomianów przepisać w formie:

$p_{1}(t)=1-{{2} \over {\beta +\alpha }}t\;$

(5.143)

wtedy macierz M możemy uzyskać podstawiając za "t" macierz A w wyrażeniu (5.143), a macierz N przepisujemy po uzyskanym wyrażeniu p₁(t) w formie:

$\mathbf {N} ={{2} \over {\beta +\alpha }}\mathbf {I} \;$

(5.144)

Wielkość ρ(M) możemy w taki sposób przepisać pamiętając przy tym, że |p(t)|<1, która jest napisane dla t∈<α,β>, by jego ograniczenie od góry przepisać w formie:

$\rho (\mathbf {M} )=1-{{2} \over {\beta +\alpha }}t\leq 1-{{2} \over {\beta +\alpha }}\alpha ={{\beta +\alpha -2\alpha } \over {\beta +\alpha }}={{\beta -\alpha } \over {\beta +\alpha }}\;$

(5.145)

Wyrażenie (5.142) przepiszemy dla s=2, przy którym dla tego s napiszemy odpowiednio wielomiany Czybyszewa (1.27), którego wielomian jest napisany w postaci rozwiniętej T₂(x)=2x²-1, które wykorzystamy do wspomnianego wyrażenia, by potem przepisać go w formie:

$p_{2}(t)={{T_{2}\left({{2t-(\beta +\alpha )} \over {\beta -\alpha }}\right)} \over {T_{2}\left(-{{\beta +\alpha } \over {\beta -\alpha }}\right)}}={{2\left({{2t-(\beta +\alpha )} \over {\beta -\alpha }}\right)^{2}-1} \over {2\left(-{{\beta +\alpha } \over {\beta -\alpha }}\right)^{2}-1}}={{2(2t-(\beta +\alpha ))^{2}-(\beta -\alpha )^{2}} \over {2(\beta +\alpha )^{2}-(\beta +\alpha )^{2}}}=\;$
$={{8t^{2}+2(\beta +\alpha )^{2}-8t(\beta +\alpha )-(\beta -\alpha )^{2}} \over {2\beta ^{2}+\alpha ^{2}+4\alpha \beta -\beta ^{2}-\alpha ^{2}+2\alpha \beta }}={{8t^{2}-8t(\beta +\alpha )+\beta ^{2}+\alpha ^{2}+6\beta \alpha } \over {\beta ^{2}+\alpha ^{2}+6\beta \alpha }}=\;$
$=1-t{{8(\beta +\alpha )} \over {\beta ^{2}+6\beta \alpha +\alpha ^{2}}}+t^{2}{{8} \over {\beta ^{2}+6\beta \alpha +\alpha ^{2}}}$

(5.146)

Wielkość ρ(M) możemy w taki sposób przepisać pamiętając przy tym, że |p(t)|<1 jest napisane dla t∈<α,β>, by jego ograniczenie od góry przepisać w formie:

$\rho (\mathbf {M} )=1-t{{8(\beta +\alpha )} \over {\beta ^{2}+6\beta \alpha +\alpha ^{2}}}+t^{2}{{8} \over {\beta ^{2}+6\beta \alpha +\alpha ^{2}}}\leq \;$
$\leq 1-\alpha {{8(\beta +\alpha )} \over {\beta ^{2}+6\beta \alpha +\alpha ^{2}}}+\beta ^{2}{{8} \over {\beta ^{2}+6\beta \alpha +\alpha ^{2}}}={{\beta ^{2}+6\beta \alpha +\alpha ^{2}-8\alpha (\beta +\alpha )+8\beta ^{2}} \over {(\beta +\alpha )^{2}+4\beta \alpha }}=\;$
$={{9\beta ^{2}-7\alpha ^{2}-2\alpha \beta } \over {(\beta +\alpha )^{2}+4\beta \alpha }}={{(\beta -\alpha )^{2}+8(\beta ^{2}-\alpha ^{2})} \over {(\beta +\alpha )^{2}+4\beta \alpha }}$

(5.147)

Patrząc na wzory (5.146) dla s=1 i na (5.147) dla s=2 dochodzimy do wniosku, że dla s=1 przy y⁽⁰⁾=x⁽ⁱ⁾ mamy:

$\mathbf {y} ^{(1)}=\left(\mathbf {I} -{{2} \over {\beta +\alpha }}{\vec {A}}\right)\mathbf {x} ^{(i)}+\mathbf {b} =\mathbf {x} ^{(i)}-{{2} \over {\beta +\alpha }}\left(\mathbf {A} \mathbf {x} ^{(i)}-\mathbf {b} \right)\;$

(5.148)

a także dla dowolnego k można powiedzieć dla k=1,2,..,s-1:

$\mathbf {y} ^{(k+1)}=\mathbf {y} ^{(k)}+\omega _{k}\omega _{k-1}(\mathbf {y} ^{(k)}-\mathbf {y} ^{(k-1)})-{{2} \over {\beta +\alpha }}(1+\omega _{k}\omega _{k-1})-(\mathbf {A} \mathbf {y} ^{(k)}-\mathbf {b} )\;$

(5.149)

We wzorach (5.148) i (5.149) możemy zdefiniować liczby ω₀, ω₁,... , które wyliczamy jako:

$\omega _{0}={{\beta -\alpha } \over {\beta +\alpha }}\;$

(5.150)

$\omega _{k+1}={{1} \over {2{{\beta +\alpha } \over {\beta -\alpha }}-\omega _{k}}}\;$

(5.151)

Twierdzenie
Określmy macierz dodatnio określoną i dodatkowo symetryczną, dla której wartości t są określone dla przedziału <α,β> przy którym α jest wartością większą niż zero, tak by dla dowolnej wartości początkowej wektora x₀ i wektora x⁽¹⁾,x⁽²⁾,..., które są generowane przez ogólny wzór według algorytmu (5.149), który jest zbieżny do ogólnego rozwiązania Ax=b.

Dowód
Algorytm (5.149) wynika z (5.142), która wynika z p(t)=1-w(t)t, który w miejsce t podstawiamy A, wtedy powstaje M, to ciąg (5.133) generuje ciąg wektorów zbieżnych do dokładnego rozwiązania Ax=b.