Maticový počet – základní pojmy Matice je obdélníkové schéma tvaru a11
a12 . . . ...... a1n
a21
a22 …………. a2n , kde aij ∈ R ( i = 1,…,m, j = 1,…,n )
A= am1 • • • • • • • • •
am2 ………. amn
aij ∈ R se nazývají prvky matice o matici o m řádcích a n sloupcích říkáme, že je typu m/n a značíme A(m/n) je-li aij = 0 pro ∀ i = 1…m, j = 1…n nazýváme matici nulovou je-li počet řádků roven počtu sloupců (tj. m=n) nazýváme příslušnou matici A(n/n) čtvercovou maticí řádu n. prvky matice aij, pro které platí i=j, tvoří hlavní diagonálu matice (platí pro čtvercovou matici) prvky matice aij, pro které platí j=n-i+1, tvoří vedlejší diagonálu matice (platí pro čtvercovou matici řádu n) čtvercová matice řádu n, která má na hlavní diagonále samé 1 a ostatní prvky jsou nulové se nazývá jednotková matice čtvercová matice řádu n, která má mimo hlavní diagonálu všechny prvky nulové se nazývá diagonální matice matice A’ (n / m), která vznikne z matice A (m / n) záměnou řádků za sloupce se nazývá transponovaná matice
Numerické řešení soustavy lineárních rovnic Řešení soustav lineárních rovnic patří mezi nejdůležitější části numerické matematiky. Mnoho praktických úloh nakonec vede k řešení takovýchto soustav, často velmi rozsáhlých. Budeme se zabývat řešením soustavy n lineárních rovnic a11 x1 a21 x1 : an1 x1
+ a12 x2 + … + a1n xn = a1,n+1 + a22 x2 + … + a2n xn = a2,n+1 : + an2 x2 + … + ann xn = an,n+1
(1)
s n neznámými x1, x2, ..., xn. Soustavu můžeme psát ve tvaru
ai,n+1 nebo v maticovém tvaru A.x=b
(2)
kde
a1,n+1 a2,n+1 :
an,n+1 Matice A=(aij), i,j=1,..., n , se nazývá matice soustavy a sloupcový vektor b=(a1,n+1,...,an,n+1) vektor pravých stran a x=( x1,...,xn) je vektor neznámých.
a11 a12 ……… a1n a21 a22 ……… a2n . . an1 an2 ……… ann
matice soustavy
a11 a12 ……… a1n a21 a22 ……… a2n . . an1 an2 ……… ann
a1,n+1 a2,n+1
an,n+1
rozšířená matice soustavy
Řešením soustavy rozumíme uspořádanou n-tici (r1, r2, …, rn), kde po dosazení rj za xj do všech rovnic přejdou tyto rovnice v identitu. Soustava lineárních rovnic se nazývá řešitelná (resp. neřešitelná), jestliže existuje (neexistuje) alespoň jedno její řešení. Dvě soustavy nazýváme ekvivalentní, jestliže množiny jejich řešení jsou si rovny. Jakoukoliv úpravu soustavy lineárních rovnic, po níž vznikne soustava ekvivalentní, nazýváme ekvivalentní úpravou. Ekvivalentní úpravy: • libovolná záměna pořadí rovnic • vynásobení libovolné rovnice nenulovým číslem ze stejné množiny jako jsou koeficienty • přičtení libovolné lineární kombinace ostatních rovnic k rovnici • přidání rovnice, která je lineární kombinací rovnic soustavy • vypuštění rovnice, která je lineární kombinací zbývajících rovnic Je-li matice soustavy regulární má soustava právě jedno řešení v R, je-li matice soustavy singulární, má soustava v R nekonečně mnoho řešení nebo nemá žádné. Všude v dalším textu budeme předpokládat, že matice soustavy je regulární, tj. že řešená soustava má právě jedno řešení. Metody pro řešení soustav lineárních rovnic dělíme na přímé a iterační.
Přímé metody Přímé metody vedou k řešení soustavy po konečném počtu kroků. Takto nalezené řešení by bylo přesné, kdybychom se v průběhu výpočtu nedopouštěli zaokrouhlovacích chyb. 1. Gaussova eliminační metoda (GEM) Základem této metody je úprava soustavy na trojúhelníkový tvar pomocí ekvivaletních úprav. Tuto část GEM nazýváme přímý chod GEM. Přidáme-li v soustavě (2) vektor pravých stran b jako (n+1)-ní sloupec k matici A, můžeme soustavu přepsat ve tvaru a11 x1 + a12 x2 + .. . + a1n xn = a1,n+1 a21 x1 + a22 x2 + .. . + a2n xn = a2,n+1 : : an1 x1 + an2 x2 + .. . + ann xn = an,n+1 Nyní se pomocí přičítání vhodných násobků první rovnice budeme snažit z ostatních rovnic eliminovat x1. Je-li a11=0, vyměníme první rovnici s první takovou rovnicí, která na prvním místě nulu nemá. Odečteme-li postupně první rovnici, vynásobenou číslem [(ai1)/(a11)], od i-té rovnice, pro i=2,3, ... , n, dostaneme a11 x1
+ a12 x2 + .. . + a1n xn = a1,n+1 a22 x2 + .. . + a2n xn = a2,n+1 : : an2 x2 + .. . + ann xn = an,n+1
Nové koeficienty jsou vypočteny jako aij=aij - [(ai1)/(a11)] . a1j , i = 2,3,..., n,
j=2,3,..., n+1.
Nyní budeme pomocí vhodných násobků druhé rovnice eliminovat x2 v třetí, čtvrté, ... n-té rovnici. Opět, je-li a22=0, vyměníme druhou rovnici s první z dalších rovnic, ve které u x2 nula není. Tím dostaneme a11 x1 + a12 x2 + a13 x3 +.. . + a1n xn = a1,n+1 a22 x2 + a23 x3 +.. . + a2n xn = a2,n+1 a33 x3 +.. . + a3n xn = a3,n+1 : : an3 x3 +.. . + ann xn = an,n+1
kde
aij=aij -[(ai2)/(a22)] . a2j ,
i = 3,4,..., n,
j=3,4,..., n+1.
Pokračujeme-li dále stejným způsobem, dostaneme po n-1 krocích soustavu v trojúhelníkovém tvaru a11 x1 + a12 x2 + a13 x3 +.. . + a1n xn = a1,n+1 a22 x2 + a23 x3 +.. . + a2n xn = a2,n+1 a33 x3 +.. . + a3n xn = a3,n+1 : : ann xn = an,n+1
(3)
Z této soustavy snadno určíme hledané řešení:
x x
n
a a
n , n +1
=
n −1
x= 1
=
n,n
(
1
a
1
a
n −1, n −1
. a n −1,n +1 − a n−1, n xn
)
(4)
(
. a1,n +1 − a12 x2 − a13 x3 − ... − a1n xn
)
11
Postup vedoucí k soustavě (3) se nazývá přímý chod GEM, výpočet neznámých dle (4) zpětná substituce nebo též zpětný chod GEM. Číslo akk nazýváme hlavní prvek. Příklad: Pomocí Gaussovy eliminace vyřešte soustavu rovnic 1,67 x - 0,15 x + 2,51 x = -0,84 2,15 x + 3,02 x - 0,17 x = 2,32 1,71 x - 2,83 x + 1,45 x = 1,26 1
2
1 1
3
2
3
2
3
Řešení: Koeficienty soustavy opíšeme do matice: 1,6 7 2,1 5 1,7 1
-0, 15 2,5 1 -0, 84 3,0 2 -0, 17 2,3 2 -2, 83 1,4 5 1,2 6
Od druhého řádku odečteme první řádek vynásobený [2,15/1,67] a od třetího vynásobený [1,71/1,67] (všechny mezivýsledky jsou zaokrouhlovány na pět desetinných míst): 1,67 0 0
-0,15 3,21311 -2,67641
2,5 1 -3,40 144 -1,12 012
-0, 84 3,40 144 2,12 012
Nyní od třetího řádku odečteme druhý vynásobený [2,67641)/3,21311]. Tím dostaneme 1,6 7 0 0
-0, 15 3,21 311 0
2,5 1 -3,40 144 -3 ,95 339
-0, 84 3,40 144 4,95 339
což už odpovídá soustavě v trojúhelníkovém tvaru 1,67 x
1
- 0,1 5 x 3,2131 1 x 2
+ 2,51 x - 3,4014 4 x - 3,9533 9 x 3
2
3 3
= - 0,84 = 3,40144 = 4,95339
Řešení této soustavy je
4,95339 ≅ −1,25295 − 3,95339 1 x2 = 3,21311 (3,40144 + 3,40144.(− 1,25295)) ≅ −0,26777
x
3
=
x
1
=
1 (− 0,84 + 0,15.(− 0,26777 ) − 2,51.(− 1,25295)) ≅ −1,35613 1,67
Řešení získané Gaussovou eliminační metodou by bylo přesné, kdybychom se v průběhu výpočtu nedopouštěli zaokrouhlovacích chyb. U některých soustav může být bohužel vliv zaokrouhlování na výsledek značný. Algoritmus Gaussovy eliminace se proto někdy modifikuje způsobem popsaným v následující kapitole.
1.1 Eliminace s výběrem hlavního prvku Eliminace s výběrem hlavního prvku je modifikace Gaussovy eliminační metody, která slouží ke zmenšení zaokrouhlovacích chyb. Abychom se vyhnuli dělení čísly, která jsou malá vzhledem k ostatním veličinám, použijeme postup zvaný výběr hlavního prvku: V prvním kroku eliminace najdeme rovnici, která má u x1 v absolutní hodnotě největší koeficient. Vyměníme ji s první rovnicí a pak pomocí jejích násobků eliminujeme x1 z ostatních rovnic. Ve druhém kroku najdeme mezi všemi rovnicemi kromě první tu rovnici, která má v absolutní hodnotě největší koeficient u x2. Vyměníme ji s druhou rovnicí a pomocí jejích násobků eliminujeme x2 z dalších rovnic. Obecně v k-tém kroku eliminace najdeme mezi posledními n-k+1 rovnicemi tu, která má největší koeficient u xk, vyměníme ji s k-tou rovnicí a pak pomocí ní eliminujeme. Právě popsanou metodu bychom mohli nazvat výstižněji eliminační metodou s částečným výběrem hlavního prvku.
Úplný výběr hlavního prvku spočívá v tom, že v k-tém kroku volíme za hlavní prvek ten, který je největší v absolutní hodnotě v submatici vytvořené vynecháním prvních k-1 řádků a sloupců v upravované matici. Nutnost hledat největší prvek v celé submatici a vyměňovat řádky i sloupce způsobuje větší časovou (a programátorskou) náročnost této metody. Gaussova eliminační metoda s částečným výběrem je proto obvykle efektivnější než metoda s úplným výběrem hlavního prvku.
Algoritmus pro přímý chod GEM (pouze část, která nuluje prvky pod hlavní diagonálou v i-tém sloupci za předpokladu, že prvek na hlavní diagonále v daném sloupci je nenulový)
j := i + 1 p := aji / aii k := 1 ajk := ajk - p * aik k := k + 1 -
k > n+1 +
j := j + 1 -
j>n +
nulování prvků v i-tém sloupci probíhá od (i+1). rovnice
Algoritmus pro zpětný chod GEM (za předpokladu, že matice soustavy je regulární)
xn := an,n+1 / an,n i := n-1 s := 0 j := i + 1 s := s + xj * aij j := j + 1 -
j>n +
xi := ai,n+1 - s / ai,i i := i - 1 -
i<1 +
2. Jordanova metoda Jordanova metoda je modifikací GEM. Přímý chod zůstává stejný, tj. ekvivalentními úpravami převedeme matici soustavy na trojúhelníkový tvar. Dalšími ekvivalentními úpravami převedeme matici soustavy na matici diagonální, tj. prvky pod i nad hlavní diagonálou jsou nulové. a11 x1 a22 x2 a33 x3 : ann xn Nejprve eliminujeme x2 u 1. rovnice aij=aij - [(ai2)/(a22)]. a2j , i = 1,
= a1,n+1 = a2,n+1 = a3,n+1 : = an,n+1
j=2,3,..., n+1.
a získáme soustavu ve tvaru a11 x1 + a13 x3 +.. . + a1n xn = a1,n+1 a22 x2 + a23 x3 +.. . + a2n xn = a2,n+1 a33 x3 +.. . + a3n xn = a3,n+1 : : ann xn = an,n+1 Dále pak eliminujeme x3 u 1. a 2. rovnice aij=aij - [(ai3)/(a33)]. a3j , i = 1,2 a získáme soustavu ve tvaru a11 x1 a22 x2
j=3,..., n+1.
.. . + a1n xn = a1,n+1 .. . + a2n xn = a2,n+1 a33 x3 + .. . + a3n xn = a3,n+1 : : ann xn = an,n+1
Soustava je nakonec ve tvaru: a11 x1 a22 x2 a33 x3 : ann xn
= a1,n+1 = a2,n+1 = a3,n+1 : = an,n+1
Řešení soustavy pak dostaneme: x1 = a1,n+1 / a11 x2 = a2,n+1 / a22
: xn = an,n+1 / ann Tedy
xi = ai,n+1 / aii,
pro i = 1, 2, …, n
Shrnutí •
GEM a Jordanova metoda vedou přímo k řešení soustavy. Kdybychom se nedopouštěli zaokrouhlovacích chyb, našli bychom pomocí těchto metod přesné řešení.
•
Základem Gaussovy eliminační metody je úprava matice soustavy na trojúhelníkový tvar. Ten dostaneme pomocí ekvivalentních úprav (přičítání vhodných násobků vybraných řádků matice k ostatním řádkům)
•
Vliv zaokrouhlovacích chyb u přímých metod může být značný, zvlášť u některých typů matic. Proto se používá tzv. eliminace s výběrem hlavního prvku.
•
Eliminační metoda je velmi náročná z časového i paměťového hlediska. Nejlépe se hodí pro nepříliš rozsáhlé soustavy s plnou maticí.
Iterační metody Mnoho praktických problémů vyžaduje řešení rozsáhlých soustav lineárních rovnic Ax=b, v nichž matice A je řídká, tj. má relativně málo nenulových prvků. Standardní eliminační metody nejsou pro řešení takových soustav vhodné, neboť v průběhu eliminace dochází postupně k zaplňování původně nenulových pozic v matici soustavy, což vede k velkým nárokům na počet aritmetických operací a klade také vysoké nároky na paměť počítače. To je důvod, proč se pro řešení takových soustav používají iterační metody. Zvolí se počáteční vektor x0 a generuje se posloupnost vektorů x0 ------> x1-----> x2 -----> … , která konverguje k hledanému řešení x. Iterační metody, na rozdíl od přímých metod, nevedou k přesnému řešení po konečném, předem daném počtu kroků. U iteračních metod zvolíme počáteční aproximaci řešení a určitým postupem ji v každém kroku metody zlepšíme. K řešení se přibližujeme postupně a obecně ho dosáhneme až v limitě. Protože výpočet nelze provádět do nekonečna, po jisté době jej ukončíme. Výsledkem bude přibližné řešení soustavy. 1. Jacobiho metoda Popíšeme, jak se Jacobiho metodou soustavy rovnic řeší a kdy se touto metodou řešit mohou. Budeme opět pracovat se soustavou lineárních rovnic a11 x1 + a12 x2 + .. . + a1n xn = a1,n+1 a21 x1 + a22 x2 + .. . + a2n xn = a2,n+1 : : an1 x1 + an2 x2 + .. . + ann xn = an,n+1 Z první rovnice vyjádříme x1, ze druhé rovnice x2 atd. Dostaneme x1 = 1/ a11 (a1,n+1 - a12 x2 - a13 x3 - .. . - a1n xn ) x2 = 1/ a22 (a2,n+1 - a21 x1 – a23 x3 -.. . - a2n xn ) : : xn = 1/ ann (an,n+1 - an1 x1 – an2 x2 -.. . - an,n-1 xn-1 )
(6)
Řešení soustavy budeme hledat následujícím způsobem: Libovolně zvolíme počáteční aproximaci řešení x(0)=(x1(0),x2(0),...,xn(0)). Tato čísla dosadíme do pravé strany (6). Tím dostaneme novou aproximaci řešení x(1)=(x1(1),x2(1),...,xn(1)). Tu opět dosadíme do pravé strany (6) atd. Obecně každou další aproximaci řešení získáme podle předpisu x1 (r+1) = 1/ a11 (a1,n+1 - a12 x2(r)- a13 x3(r)- .. . - a1n xn(r) ) x2 (r+1) = 1/ a22 (a2,n+1 - a21 x1(r)– a23 x3(r)-.. . - a2n xn(r)) : : xn (r+1) = 1/ ann (an,n+1 - an1 x1(r)– an2 x2(r)-.. . - an,n-1 xn-1(r) )
(6)
Maticově můžeme tento výpočet zapsat takto: x1 x2 : : xn
r+1
0 -a21/a22 = -an1/ann
-a12/a11 0 : : -an2/ann
-a13/a11 -a23/a22
… …
. -an3/ann
…
iterační matice nebo také
x
( r +1)
i
n
= ∑− j =1
a* x a ij ii
(r )
j
+
a a
i , n +1 ii
-a1n/a11 -a2n/a22 : : 0
*
x1 x2 : : xn
r
+
a1,n+1/a11 a2,n+1/a22 : : an,n+1/ann
Za jistých (dále popsaných podmínek) se tímto postupem budeme přibližovat k přesnému řešení soustavy. Ve výpočtu pokračujeme, dokud se nedosáhne určité předem dané přesnosti (ε), např. dokud se aproximace řešení neustálí na požadovaném počtu desetinných míst |xi(r+1) – xi(r)| < ε , pro každé i=1..n nebo dokud není překročen předem daný maximální počet kroků. Jacobiho metodou nemusíme řešení soustavy najít vždy. V některých případech posloupnost postupných aproximací k řešení soustavy nekonverguje. Uvedeme nyní podmínky, které zaručí, že metoda konverguje (tj. najdeme pomocí ní přibližné řešení). Definice: Matice A se nazývá řádkově ostře diagonálně dominantní právě tehdy, když
(neboli když je v každém řádku matice absolutní hodnota prvku na diagonále větší než součet absolutních hodnot všech ostatních prvků v onom řádku) a sloupcově ostře diagonálně dominantní právě tehdy, když
(neboli když je v každém sloupci matice absolutní hodnota prvku na diagonále větší než součet absolutních hodnot všech ostatních prvků v onom sloupci). Je-li matice soustavy (2) ostře řádkově nebo sloupcově diagonálně dominantní, Jacobiho metoda konverguje. Jestliže matice soustavy (2) není diagonálně dominantní, Jacobiho metoda konvergovat může a nemusí. Existuje podmínka pro konvergenci Jacobiho metody nutná a dostatečná (tj. pokud je splněna, metoda konverguje a pokud není splněna, metoda diverguje), jenže je pro velké matice prakticky neověřitelná. Proto, nejsme-li si jisti konvergencí metody, je vhodné stanovit maximální počet kroků a je-li překročen, výpočet ukončit s tím, že metoda diverguje. Pak je potřeba zvolit jinou metodu nebo soustavu nějak upravit. Příklad: Jacobiho metodou řešte soustavu
Řešení: Matice soustavy je diagonálně dominantní, protože platí
Proto je konvergence metody zaručena. Vypíšeme iterační vztahy:
Jako počáteční aproximaci zvolíme x=(0,0,0) Postupně získávané aproximace řešení budeme zapisovat do tabulky:
Je vidět, že posloupnost postupných aproximací konverguje k řešení soustavy (2,-2,-1). Kdybychom chtěli získat řešení s přesností ε = 0,01, mohli bychom nyní výpočet zastavit, protože
zatímco kdybychom požadovali přesnost ε = 0,001, museli bychom ve výpočtu pokračovat, protože např.
2. Gauss-Seidelova metoda Gauss-Seidelova metoda je velmi podobná metodě Jacobiho. Liší se od ní pouze v tom, že při výpočtu další aproximace řešení použijeme vždy nejnovější přibližné hodnoty x1, x2, ... , xn, které máme k dispozici. Podrobněji: x1(r+1) vypočteme stejně jako u Jacobiho metody a při výpočtu x2(r+1) je ihned použijeme (zatímco u Jacobiho metody jsme použili staré x1(r)). Při výpočtu x3(r+1) použijeme nové x1(r+1) a x2(r+1) atd. Obecně iterační vztahy vypadají takto: x1 (r+1) = 1/ a11 (a1,n+1 - a12 x2(r)- a13 x3(r)- .. . - a1n xn(r) ) x2 (r+1) = 1/ a22 (a2,n+1 - a21 x1(r+1)– a23 x3(r)-.. . - a2n xn(r)) : : xn (r+1) = 1/ ann (an,n+1 - an1 x1(r+1)– an2 x2(r+1)-.. . - an,n-1 xn-1(r+1) ) nebo také
x
( r +1)
i
i −1
= ∑− j =1
a* x a ij ii
( r +1)
j
+
n
∑−
j = i +1
a* x a ij ii
(r )
j
+
a a
i , n +1 ii
Dá se dokázat, že je-li matice soustavy (3) ostře řádkově diagonálně dominantní, Gauss-Seidelova metoda konverguje. Příklad: Gauss-Seidelovou metodou řešte tutéž soustavu, t.j.
Řešení: Již jsme ověřili, že podmínka konvergence je splněna. Vypíšeme iterační vztahy:
Jako počáteční aproximaci zvolíme opět x=(0,0,0).
Vidíme, že se k řešení soustavy přibližujeme rychleji než pomocí Jacobiho metody. I obecně se dá říci, že GaussSeidelova metoda obvykle konverguje rychleji než metoda Jacobiho. Proto se používá častěji. Další její výhodou oproti Jacobiho metodě je, že pro uložení přibližného řešení v paměti počítače nám stačí jediné pole, jehož složky postupně přepisujeme, zatímco u Jacobiho metody si musíme pamatovat pole dvě: starou a novou aproximaci řešení. Shrnutí •
Pomocí iteračních metod obvykle najdeme pouze přibližné řešení soustavy (pokud nenastane dosti nepravděpodobný případ, kdy se v některém kroku trefíme přímo do řešení).
•
Na začátku zvolíme počáteční aproximaci řešení a tu pak opakovaným dosazováním do iteračních vztahů zpřesňujeme. S výpočtem skončíme obvykle tehdy, je-li norma rozdílu po sobě jdoucích aproximací dostatečně malá.
•
Iterační metody mohou divergovat (řešení pomocí nich nemusíme najít). Zda bude metoda konvergovat, či nikoli, závisí na vlastnostech matice soustavy.
•
Iterační metody jsou vhodné pro řešení velkých soustav s řídkou maticí koeficientů. Pro řešení malého počtu rovnic vhodné nejsou, tam lépe poslouží eliminace. Řídká matice soustavy Řekneme, že matice je řídká, když počet jejích nenulových prvků je výrazně menší než počet jejích prvků. Řídké matice jsou v paměti počítače účelně reprezentovány jen pomocí svých nenulových koeficientů. Jedním z hlavních cílů efektivních algoritmů pro řešení soustav s řídkými maticemi je provádět takové kroky, aby vznikalo co nejméně nových nenulových koeficientů.
Algoritmus vyjádření iterační matice
i := 1 j := 1 aij := - aij / aii j := j + 1 -
j>n +
ai,n+1 := ai,n+1 / aii aii := 0 i := i + 1 -
i>n +
Algoritmus výpočtu i-tého přiblížení u Jacobiho metody (XP … předchozí přiblížení, XN … nově vypočítávané přiblížení)
nastavení počátečního přiblížení –> xpi = 0 (pro i=1..n)
i := 1 xni := ai,n+1 j := 1 xni := xni + xpj * aij j := j + 1 -
j>n +
i := i + 1 -
i>n +
Algoritmus výpočtu i-tého přiblížení u Gauss-Seidelovy metody (X … přiblížení)
nastavení počátečního přiblížení –> xi = 0 (pro i=1..n)
i := 1 xi := ai,n+1 j := 1 xi := xi + xj * aij j := j + 1 -
j>n +
i := i + 1 -
i>n +
