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} \;$