UNIVERZITA PARDUBICE Fakulta elektrotechniky a informatiky
Iterační metody řešení soustav lineárních rovnic Václav Baťka
Bakalářská práce 2012
Prohlášení autora Prohlašuji, že jsem tuto práci vypracoval samostatně. Veškeré literární prameny a informace, které jsem v práci využil, jsou uvedeny v seznamu použité literatury. Byl jsem seznámen s tím, že se na moji práci vztahují práva a povinnosti vyplývající ze zákona č. 121/2000 Sb., autorský zákon, zejména se skutečností, že Univerzita Pardubice má právo na uzavření licenční smlouvy o užití této práce jako školního díla podle § 60 odst. 1 autorského zákona, a s tím, že pokud dojde k užití této práce mnou nebo bude poskytnuta licence o užití jinému subjektu, je Univerzita Pardubice oprávněna ode mne požadovat přiměřený příspěvek na úhradu nákladů, které na vytvoření díla vynaložila, a to podle okolností až do jejich skutečné výše. Souhlasím s prezenčním zpřístupněním své práce v Univerzitní knihovně.
V Pardubicích dne 9. 5. 2012
Václav Baťka
Poděkování V první řadě by chtěl autor poděkovat své rodině, která ho po dobu studia podporovala, jak finančně, tak psychicky. Dále by chtěl poděkovat vedoucí Mgr. Janě Heckenbergerové Ph.D., která autora v průběhu práce směrovala ke zdárnému konci. V neposlední řadě přátelům, kteří autora po dobu studia povzbuzovali.
Anotace Práce se zabývá řešením soustav lineárních rovnic, jak přímými tak hlavně iteračními metodami. Srovnává navrženou aplikaci s Matlabem. Klíčová slova Soustava lineárních rovnic, přímé metody, iterační metody, SLR,
Title Iterative methods for solving systems of linear equations
Annotation The work deals with the solving systems of linear equations, both direkt and iterative methodes mainly. Compares the proposed application with Matlab Keywords systems of linear equations, direct methods, iterative methods,
Obsah Seznam grafů ....................................................................................................................... 8 Seznam tabulek .................................................................................................................... 9 Seznam obrázků ................................................................................................................. 10 Úvod .................................................................................................................................... 11 1
Základní pojmy z teorie systémů rovnic ................................................................ 12 1.1 Soustava lineárních rovnic........................................................................................ 12
2
Přímé metody řešení systémů lineárních rovnic ..................................................... 15 2.1 Gaussova eliminační metoda .................................................................................... 15 2.2 Gaussova eliminační metoda s částečným výběrem hlavního prvku ....................... 17 2.3 Gaussova eliminační metoda s úplným výběrem hlavního prvku ............................ 19 2.4 LU rozklad ................................................................................................................ 20 2.5 Choleského metoda .................................................................................................. 22 2.6 Croutova metoda ...................................................................................................... 24
3
Nepřímé metody řešení systémů lineárních rovnic ................................................ 27 3.1 Princip iteračních metod ........................................................................................... 27 3.2 Jacobiova iterační metoda ........................................................................................ 27 3.3 Gaussova-Seidelova iterační metoda ........................................................................ 30 3.4 Ukončovací kritéria iteračních metod ...................................................................... 32
4
Algoritmy pro výpočet iteračních metod v Matlabu .............................................. 33 4.1 BICG – Metoda bi-sdružených gradientů ................................................................. 33 4.2 BICGSTAB – Stabilizovaná metoda bi-sdružených gradientů ................................ 34 4.3 CGS – Čtvercová metoda sdružených gradientů ...................................................... 34 4.4 GMRES – Zobecněná metoda minimálního rezidua ................................................ 35 4.5 QMR – Metoda částečně minimálního rezidua ........................................................ 35
5
Popis praktické části bakalářské práce ................................................................... 36 5.1 Popis grafického rozhraní programu ........................................................................ 36 5.2 Popis použitých algoritmů ........................................................................................ 37
6
Srovnání s matematickým softwarem Matlab ........................................................ 39 6.1 Srovnání časových nároků při výpočtu .................................................................... 39 6.2 Srovnání přesnosti výsledků ..................................................................................... 44
Závěr ................................................................................................................................... 46
Literatura ........................................................................................................................... 47 Příloha A – Dokumentace ................................................................................................. 49 Příloha B – Zdrojové kódy................................................................................................ 50
Seznam grafů Graf 1 Srovnání časové náročnosti Jacobiovy metody ........................................................ 41 Graf 2 Srovnání časové náročnosti Jacobiovy metody maticovým způsobem ................... 42 Graf 3 Srovnání časové náročnosti G-S metody ................................................................. 43 Graf 4 Srovnání časové náročnosti G-S metody maticovým způsobem ............................. 44
8
Seznam tabulek Tabulka 1– Iterační výpočet SLR pomocí Jacobiovy metody............................................. 30 Tabulka 2– Iterační výpočet SLR pomocí G-S metody ...................................................... 31 Tabulka 3 Srovnaní časové náročnosti Jacobiovy metody .................................................. 40 Tabulka 4 Srovnaní časové náročnosti Jacobiovy metody maticovým způsobem .............. 41 Tabulka 5 Srovnaní časové náročnosti G-S metody ............................................................ 42 Tabulka 6 Srovnaní časové náročnosti G-S metody maticovým způsobem ....................... 43 Tabulka 7 Srovnání přesnosti výsledků Jacobiovy metody................................................. 44 Tabulka 8 Srovnání přesnosti výsledků G-S metody .......................................................... 45
9
Seznam obrázků Obrázek 1 Struktura praktické části BP ............................................................................... 36 Obrázek 2 Grafické rozhraní BP ......................................................................................... 36 Obrázek 3 Kód pro výpočet Jacobiovy metody v Matlabu ................................................. 39 Obrázek 4 Kód pro výpočet G-S metody v Matlabu ........................................................... 40
10
Úvod Tato bakalářská práce si klade za cíl seznámit čtenáře s hlavními aspekty teorie systémů lineárních rovnic. V první kapitole seznamuje čtenáře se základními pojmy teorie systému lineárních rovnic. V druhé kapitole se nejdříve obecně věnuje problematice přímých iteračních metod. Poté se důkladně projdou jednotlivé přímé metody (Gaussova eliminační metoda bez výběru i s výběrem hlavního prvku, LU rozklad, Choleshéko metoda, Croutova metoda pro třídiagonální matice) a ukáže se jejich postup řešení na příkladech. V třetí kapitole se seznámí čtenář s iteračními metodami řešení soustav lineárních rovnic. Projdou se Jacobiova a Gauss-Seidelova metoda, ukáže se názorně průběh řešení na vzorových příkladech. Dále se čtenář seznámí s ukončujícími kritérii pro iterační výpočty. Čtvrtou kapitolou si čtenář projde algoritmy pro řešení soustav lineárních rovnic iteračními metodami, které jsou součástí jádra Matlabu. Pátá kapitola seznamuje čtenáře s praktickou částí bakalářské práce, ukazuje návrh grafického rozhraní aplikace pro řešení soustav lineárních rovnic iteračními metodami a dále popisuje algoritmy, které byly napsány v programovacím jazyce Java. Šestá kapitola srovnává časovou náročnost výpočtu praktické části bakalářské práce a Matlabu. Srovnává se na posloupnosti příkazů, které se vykonávají v praktické části bakalářské práce a dále srovnává přesnost výsledků.
11
1 Základní pojmy z teorie systémů rovnic Matice nad tělesem P je zobrazení {1, 2, …, n}x{1, 2, …, n}
. Zapisuje se většinou
velkými písmeny. A = {…}. Je to struktura tvořená n sloupci a m řádky. Čtvercová matice je matice tvořena n sloupci a n řádky.
Nulová matice je taková, kterou má na všech pozicích aij hodnotu 0.
Jednotková matice je taková čtvercová matice, která má na hlavní diagonále prvky s hodnotou 1 a na všech zbylých pozicích nulové prvky. Hlavní diagonála obsahuje prvky aij pro i=j.
Matice transponovaná k matici A se značí AT a platí pro ni aij = aTji. To znamená, že prvky které byly na i-tém řákdu a j-tém sloupci budou nově na j-tém řádků a i-tém sloupci.
1.1 Soustava lineárních rovnic Za soustavou lineárních rovnic považujeme libovolně velkou množinu lineárních rovnic. Obecně může být soustava m lineárních rovnic s n proměnnými zapsána následovně.
12
Proměnné x1 až xn jsou neznámé a aij jsou koeficienty vystupující v soustavě rovnic. Členy bi jsou absolutní členy soustavy nebo též vektorem pravé strany. V obecném případě mohou být koeficienty i absolutní členy komplexními čísly. Koeficienty soustavy lze zapsat také v maticovém tvaru. Tuto matici poté nazýváme jako matici soustavy. Stejně tak můžeme neznámé, ale také koeficienty pravé strany zapsat jako vektory. Vektor neznámých
a vektor pravé strany .
Pokud máme soustavu lineárních rovnic vyjádřenou jako matici soustavy, vektor neznámých a vektor pravé strany lze celou soustavu lineárních rovnic zapsat ve tvaru.
Matice A a vektor pravé strany se někdy zapisuje do jedné matice, kterou označujeme A|b. Tvoří ji celá matice A do které za poslední sloupec přidáme řádky vektoru pravé strany b.
Matice A|b se nazývá rozšířená matice soustavy. Je-li vektor pravé strany b nulový, potom se jedná o soustavu homogenní. Pokud je aspoň jeden prvek vektorů b i, i = 1, 2, ..., m nenulový, potom se jedná o soustavu nehomogenní.
13
Soustavu můžeme označit řešitelnou či neřešitelnou, jestli existuje alespoň jedno, respektive neexistuje žádné řešení. Aby soustava měla řešení, plně postačuje rovnost hodnosti matice soustavy h(A) a hodnost rozšířené matice soustavy h(A|b). Řešení nehomogenních soustav má tři konce. Prvním koncem je, že daná soustava nemá řešení, dále pak soustava má právě jedno řešení. Poslední eventualitou je soustava s nekonečným počtem řešení.
14
2 Přímé metody řešení systémů lineárních rovnic Přímé metody lze aplikovat na regulárních čtvercových maticích soustav řádu n(A є Mn), kde b je vektor pravé strany. Systém má jediné řešení.
Jsou to metody, které po aplikaci konečným počtu kroků dají přesný výsledek x *, ovšem za předpokladu že všechny aritmetické operace proběhly matematicky přesně.
2.1 Gaussova eliminační metoda Jedná se o jednu z nejznámějších přímých metod řešení systémů rovnic. “Hlavní myšlenka této metody spočívá v převedení daného systému A x = b vhodnými ekvivalentními úpravami na systém R x = c s horní trojúhelníkovou maticí R, tj. problém se převede na redukovaný problém.“1 Gaussovou eliminační metodou lze také řešit výpočet inverzní matice, nebo výpočet determinantu matice. Matice A = (aij) tvaru m x n je v řádkově odstupňovaném tvaru, jestliže splňuje dvě podmínky. Obsahuje-li i-tý řádek Ai* samé nuly, pak všechny pod i-tým řádkem jsou také nulové. Je-li v i-tém řádku Ai* první nenulový prvek v j-tém sloupci, tedy na místě aij, pak všechny prvku ležící v prvních j sloupcích a současně pod i-tým řádkem se rovnají nule. První nenulový prvek v každém řádku nazýváme hlavním prvkem. V prvním kroku algoritmu vynásobíme první řádek číslem –k1, které se součtem druhého řádku odstraní hodnotu a21 ze soustavy. V druhém kroku taktéž vynásobíme první řádek číslem –k2 které tentokrát eliminuje hodnotu a31. Ve třetím kroku vynásobíme třetí řádek –k3 a tím odstraníme poslední prvek a32 k dokončení algoritmu.
1
(Horová, a další, 2008)
15
Následující krok rozhoduje o řešitelnosti soustavy. Soustavy řešené Gaussovou eliminační metodou jsou neřešitelné, pokud h(A|b) > h(A), vždy když případe rovnosti hodností matic h(a|b)=h(a).
Jsou řešitelné
V kroku tři musíme zjistit počet neznámých, které si můžeme zvolit a počet neznámých vyjádřených parametrem.
K je číslo volnosti, které nám ukazuje kolik neznámých si můžeme zvolit, nebo které musíme nahradit parametrem. N je celkový počet neznámých obsažených v soustavě. V případě kdy dostaneme k=0, tak výsledek soustavy je jednoznačný - například x1 = 1, x2 = 2, x3 = 3. Při výsledku, kdy se k rovná jedné, výsledek soustavy bude v podobném tvaru x1 = 1 + p, x2 = 2 - p, x3 = p. Kde p si zvolíme a soustava má nekonečně mnoho řešení. V posledním kroku algoritmu z redukované matice A|b se vyjádříme neznámé.
Nyní si tuto metodu názorně převedeme na příkladu.
Eliminujeme hodnoty a21, a31 na hodnotu 0. První řádek násobíme a odečteme od druhého řádku. Poté první násobíme 2 a odečteme od třetího řádku.
16
Druhý řádek již máme ve tvaru, který potřebujeme. Nyní druhý řádek vynásobíme -7 a sečteme s posledním řádkem rozšířené matice soustavy.
V tuto chvíli se nám soustava natolik zjednodušila do počitatelného tvaru.
K řešení neznámých využijeme tzv. zpětný chod. Při této metodě se výpočet začíná od poslední neznámé.
2.2 Gaussova eliminační metoda s částečným výběrem hlavního prvku
Pří výpočtu Gaussovou eliminační metodou musíme předpokládat, že v průběhu algoritmu nám vznikají nenulová čísla, kterými v dalších krocích dělíme. Protože jsou tato čísla velmi malá, v absolutní hodnotě blížící se nule, vznikají při zaokrouhlování velké nepřesnosti a tím ovlivňují celkový výsledek. Toto se dá řešit metodami s výběrem hlavního prvku. V první řadě je to řádkový výběr hlavního prvku. V k-tém kroku prohledáváme k-tý řádek a hledáme prvek, který je v absolutní hodnotě větší, než zbylé hodnoty v řádku. 17
Prohledáváme sloupce, ze kterých jsme dosud nevybrali hlavní prvek. Hlavní prvek se nachází v p-tém řádku a q-tém a platí pro něj
Druhou možností je sloupcový výběr hlavního prvku. Kdy v k-tém kroku vybíráme z ktého sloupce hodnotu, který je v absolutní hodnotě maximální ze všech hodnot v daném sloupci. Výběr probíhá v řádcích, ze kterých ještě nebyl vybrán hlavní prvek. Hlavní prvek je potom v p-tém řádku a platí pro něj
Následující příklad je řešen Gaussovou eliminační metodou s částečným výběrem hlavního prvku ve sloupcích.
18
2.3 Gaussova eliminační metoda s úplným výběrem hlavního prvku
Základní myšlenkou metody je prohledávat submatici, která se zmenšuje každým již hotovým řádkem trojúhelníkové matice. Jistou nevýhodou eliminace s úplným výběrem hlavního prvku je potřeba měnit pořadí jak řádků, tak sloupců, tím je tato metoda méně efektivnější než eliminace s částečným výběrem hlavního prvku. Hlavní prvek v k-tém kroku eliminace, k=1, …, n-1, vybíráme ze sloupců a řádků, ze kterých jsme dosud nevybírali. Hlavní prvek se nachází v p-tém řádku a q-tém sloupci a platí pro něj
Nejprve prohledáváme celou matici a hledáme prvek, který je v absolutní hodnotě maximální z celé soustavy.
Nyní hledáme prvek od indexu a22 v našem případě je to posunovat nemusíme, vystačíme si pouze výměnou sloupců.
19
. Tentokrát ale již řádky
V průběhu eliminace jsme zaměňovali sloupce, proto teď na konci algoritmu musíme opět provézt zpětně výměny sloupců vektoru výsledků
2.4 LU rozklad Alan Turing, autor LU rozkladu matice postupoval tak, že postupně počítal mtý řádek matice U a m-tý sloupce matice L, m = 1, …, n. Příslušné vzorce odvozoval pomocí vzorce pro násobení matic.
„Je-li matice řádkově, nebo sloupcově diagonálně dominantní tj.
nebo
pak LU rozklad existuje. Speciálně, je-li matice sloupcově diagonálně dominantní, je “2
2
(Felcman, 2012)
20
Pro regulární matice, ve kterých lze měnit pořadí řádků, tento rozklad lze udělat. Matice U je horní trojúhelníková matice, kterou dostaneme výše zmíněnou Gaussovou eliminací. Matice L je dolní trojúhelníková a na hlavní diagonále jedničky.
Vše si názorně předvedeme na příkladu.
Nyní si naplníme hodnoty matic L a U.
21
V tuto chvíli máme naplněné matice L a U a můžeme pokračovat dále ve výpočtu. V první řadě musíme vypočítat soustavu dosadíme do další soustavy
, dále z této soustavy dostaneme vektor , který a tím dostaneme vektor výsledku
V tuto chvíli přistupujeme k poslednímu kroku.
2.5 Choleského metoda Autor metody André-Louis Cholesky ji založil na předpokladu pozitivně definitivní čtvercové matici A a na součinu horní a dolní trojúhelníkové matici a druhá matice je z první transponovaná. Matice L se nazývá Choleského trojúhelník matice A. Využití tato metoda má při výpočtu inverzní matice podle vzorce , ale hlavně se využívá při řešení soustav lineárních rovnic.
22
Metoda opět směřuje k řešení dvou systémů trojúhelníkových matic Výpočet probíhá po sloupcích, vzorce pro výpočet jednotlivých prvků vypadají následujícím způsobem.
Průběh výpočtu si ukážeme na následujícím případě, kdy máme zadanou soustavu rovnic.
Teď když už máme připravenou matici L tak smíme přistoupit k dalšímu výpočtu
23
.
Nyní získáme hodnotu vektoru (1,
, 2) a pokračujeme v řešení posledním krokem
a tím získáme hodnotu vektoru neznámých .
2.6 Croutova metoda Prescott Durand Crout navrhl metodu, která dokáže třídiagonální matice převést na součin dolní a horní trojúhelníkové matice.
Potřebujeme matici A nyní rozložit na matice L a U
Je tedy potřeba vyjádřit koeficienty matice L a také matice U.
24
Z těchto uvedených vzorců si již snadno dopočteme všechny potřebné koeficienty do obou matic. Poté co vyjádříme matici A jako A= L U, poté systém Ax = b převedeme na již známý systém
Nyní vypočítáme vektor
. Demonstrujme tuto metodu na následujícím příkladě.
ze kterého následně zjistíme řešení této soustavy tedy vektor .
25
26
3 Nepřímé metody řešení systémů lineárních rovnic Tyto metody nazýváme také metodami iteračními. Při aplikaci těchto metod nedostaneme přesná, ale aproximovaná řešení. Využívají se v případě, když máme matice velkých řádů, při kterých by výpočet přímými metodami byl velmi zdlouhavý a paměťově náročný.
3.1 Princip iteračních metod Základním rysem iteračních metod je vyjádřit soustavu lineárních rovnic ve tvaru.
„x* je řešení systému právě tehdy, když x* je řešením systému , * -1 0 n x = (E - T) g za předpokladu, že E – T je regulární. Nechť x R je libovolná počáteční aproximace. Posloupnost určená rekurentně vztahem
se nazývá iterační posloupnost a matice T se nazývá iterační matice.“3
3.2 Jacobiova iterační metoda Tvůrcem metody byl pruský matematik Carl Gustav Jacob Jacobi. Matici soustavy můžeme rozepsat do následného tvaru
kde D je matice diagonální, L je dolní trojúhelníková matice, U horní trojúhelníková matice.
3
(Horová, a další, 2008)
27
Je potřeba tuto rovnici transformovat.
Vycházíme z předpokladu, že jednotlivé prvky regulární a lze tudíž vypočítat
1, …, n, je matice D
a z tohoto vztahu získáme maticový tvar Jacobiovy iterační metody.
Lze také si také touto metodou zapsat vektor
následovně.
„Pro některé speciální matice A je zaručena konvergence Jacobiovy iterační metody. a) Silné řádkové sumační kriterium: Nechť matice A je ryze řádkově diagonálně dominantní, tj.
Pak Jacobiova iterační metoda konverguje pro každou počáteční aproximaci . b) Silné sloupcové sumační kriterium: Nechť matice A je ryze sloupcově diagonálně dominantní, tj.
Pak Jacobiova iterační metoda konverguje pro každou počáteční aproximaci .“4
4
(Horová, a další, 2008)
28
Příklad:
Nyní z rovnice převedeme na iterační tvar, tak že z každé rovnice vyjádřím jednu neznámou.
Pokud si zvolíme počáteční stav výsledek po průběhu deseti iterací je
a hodnotu konečného kritéria .
29
, potom
Tabulka 1– Iterační výpočet SLR pomocí Jacobiovy metody
3.3 Gaussova-Seidelova iterační metoda Metoda je velmi podobná předešlé Jacobiově metodě, s tím rozdílem, že pro počítá již s vyčíslenými koeficienty .
Maticový zápis se dá zapsat takto.
Za předpokladu, že
, je matice D – L regulární.
Potom je Gaussova-Seidelova iterační metoda tvaru.
30
, i > 1,
„Posloupnost generovaná Gaussovou-Seidelovou iterační metodou konverguje pro každou počáteční aproximaci právě tehdy, když “5
Metodu provedeme na stejném příkladě jako Jacobiovu, opět s počátečními podmínkami a hodnotu konečného kritéria .
Tabulka 2– Iterační výpočet SLR pomocí G-S metody
Při porovnání obou tabulek výsledků dojdeme k závěru, že výsledek z desáté iterace Jacobiovou a výsledek sedmé iterace G-S metodou se rovnají. Stačilo nám k tomu o tři iterace méně, díky použití již vypočítaných hodnot. . 5
(Horová, a další, 2008)
31
3.4 Ukončovací kritéria iteračních metod Ukončovacím kritériem je myšlen předpis, kterým získáme reálné číslo. Toto číslo následně porovnání s hodnotou tolerance odchylky dvou posledních iterací. Pokud je toto číslo větší než tolerance, pokračuje výpočet další iterací. V opačném případě je iterační proces ukončen a výsledek systému rovnic je poslední iterace. „Nechť resp. je vektorový prostor všech uspořádaných n-tic komplexních resp. reálných čísel. Prvky tohoto prostoru budeme zapisovat ve sloupcových vektorech. Vektorová norma
je funkce
(z
) s následujícími vlastnostmi:
do
1) 2) 3) 4)
“6
Norma vektoru udává velikost vektoru , kterou získáme
.
Tímto předpisem se označuje absolutní chyba, která bere rozdíl vektorů posledních dvou iterací, ze kterých zjistíme normu vektoru a tu porovnáme s , což je hodnota tolerance.
Relativní chyba se počítá velmi podobně jako absolutní chyba. Rozdíl je v tom, že ještě normu vektoru vydělíme normou vektoru .
Residuum se získá z tolerance.
6
, které musí pro ukončení iteračního výpočtu menší než
(Horová, a další, 2008)
32
4 Algoritmy pro výpočet iteračních metod v Matlabu Pro řešení iteračních metod v Matlabu je zapotřebí mít čtvercovou matici. Metody jsou navrženy pro řešení systému Ax = b, nebo pro minimalizaci výrazu || b – Ax||. Algoritmy, které Matlab používá k řešení systému lineárních rovnic iterační metodou, jsou: BICG – Biconjungate gradient method (metoda bi-sdružených gradientů) BICGSTAB – Biconjungate gradient stabilized method (stabilizovaná metoda bisdružených gradientů). CGS – Conjungate gradients square method (čtvercová metoda sdružených gradientů) GMRES – Generalized minimum residual method (zobecněná metoda minimálního rezidua) QMR – Quasi-minimal residual method (metoda částečně minimálního rezidua). Všechny metody využívají předpoklady M. Systém lineárních rovnic Ax = b je pak zaměněn za ekvivalentní systém . Předpoklad M je zvolen pro zrychlení konvergence iteračních metod. Většina iteračních metod si při řešení iteračních rovnic používá funkci \ tzv. backslash. A\B je matice podílu A a B, který je zhruba stejný jako inv(A)*B, ale je počítán jiným postupem. Jestliže A je matice typu nxn a B je sloupcový vektor o n složkách, nebo matice složená z několika takových sloupů, pak X = A\B je řešením rovnice A*X=B Gaussovou eliminací. Varovná zpráva se zobrazí, je-li A špatného typu, nebo je-li A skoro singulární. Jestliže A je matice typu mxn a B je sloupcový vektor o m složkách, nebo matice složená z takových sloupců, pak X = A\B je řešením systému A*X=B ve smyslu metody nejmenších čtverců.
4.1 BICG – Metoda bi-sdružených gradientů X = BICG(A, B) Pokouší se řešit systém lineárních rovnic A*X = B pro X. Matice koeficientů A musí být čtvercová typu nxn a sloupcový vektor musí být délky n. BICG začíná iterační proces v počátečním odhadu, který je standardně nulový vektor délky n. Iterace se opakují, dokud metoda nekonverguje, nebo dokud není spočítán maximální počet iterací. Konvergence je dosaženo, jestliže relativní reziduum NORM(B-A*X)/NORM(B) je menší nebo rovno toleranci metody. Tolerance této metody je 1e-6. Standardní počet iterací je od n do 20. Ne BICG(A, B, TOL) určuje toleranci metody. BICG(A, B, TOL, MAXIT) určuje maximální počet iterací.
33
BICG(A, B, TOL, MAXIT, M) a BICG(A, B, TOL, MAXIT, M1, M2) používají předpoklad M nebo M = M1 * M2 a efektivně řeší systém rovnic inv(M)*A*X = inv(M)*B pro X. BICG(A, B, TOL ,MAXIT ,M1 ,M2 ,X0) specifikuje počáteční odhad. Jestliže X0 je zadán jako prázdná matice, pak je standardně použit nulový vektor. X = BICG(A, B, TOL, MAXIT, M1 ,M2, X0) vrací řešení X. jestliže BICG konverguje, pak se zobrazí informační zpráva. Jestliže u BICG selhala konvergence po maximálním množství iterací nebo se BICG zastavil z jakéhokoliv důvodu, pak se zobrazí varovná zpráva informující o relativním reziduu a o čísle iterace, při kterém se metoda zastavila či selhala.
4.2 BICGSTAB – Stabilizovaná metoda bi-sdružených gradientů Funkce BICGSTAB používá stejnou syntaxi jako funkce BICG. X = BICGSTAB(A, B) BICGSTAB(A, B, TOL) BICGSTAB(A, B, TOL, MAXIT) BICGSTAB(A, B, TOL, MAXIT, M) BICGSTAB(A, B, TOL, MAXIT, M1, M2) BICGSTAB(A, B, TOL, MAXIT, M1, M2, X0) X = BICGSTAB(A ,B, TOL, MAXIT, M1, M2, X0)
4.3 CGS – Čtvercová metoda sdružených gradientů Využívá stejnou syntaxi jako BICG. X = CGS(A, B) CGS (A, B, TOL) CGS (A, B, TOL, MAXIT) CGS (A, B, TOL, MAXIT, M) CGS (A, B, TOL, MAXIT, M1, M2) CGS (A, B, TOL, MAXIT, M1, M2, X0) 34
X = CGS (A, B, TOL, MAXIT, M1, M2, X0)
4.4 GMRES – Zobecněná metoda minimálního rezidua X = GMRES(A, B) a X = GMRES(A, B, RESTART) řeší systém lineárních rovnic A*X = B pro X. A musí být čtvercová matice typu nxn a sloupcový vektor B musí být délky n. Na rozdíl od BICG, se GMRES sám restartuje při dosažení iterace RESTART, kterou použije jako nový počáteční odhad. Standardní hodnota restart je n. Konvergence je dosaženo, jestliže relativní reziduum NORM(B-A*X)/NORM(B) je menší nebo rovno toleranci metody. Tolerance této metody je 1e-6. Standardní maximální počet iterací je menší z čísel n/RESTART a 10. GMRES(A, B, RESTART, TOL) určuje toleranci metody. GMRES(A, B, RESTART, TOL ,MAXIT) určuje maximální počet iterací. GMRES(A, B, RESTART, TOL, MAXIT, M) a GMRES(A, B, RESTART, TOL, MAXIT, M1, M2) využívá předpoklad M nebo M = M1 * M2 a efektivně řeší systém rovnic inv(M)*A*X = inv(M)*b pro X. GMRES(A, B, RESTART, TOL, MAXIT, M1, M2, X0) specifikuje počáteční odhad. X = GMRES(A, B, RESTART, TOL, MAXIT, M1, M2, X0) vrací řešení X. jestliže GMRES konverguje, pak se zobrazí informační zpráva. Jestliže u GMRES selhala konvergence po maximálním množství iterací nebo se GMRES zastavil z jakéhokoliv důvodu, pak se zobrazí varovná zpráva informující o relativním reziduu a o čísle iterace, při kterém se metoda zastavila či selhala.
4.5 QMR – Metoda částečně minimálního rezidua Syntaxe příkazů je stejná jako u funkce BICG s výjimkou QMR(A, B, TOL, MATIX, M1) a QMR(A, B, TOL, MAXIT, M1,M2). U těchto příkazů je použito levého předpokladu M1 a pravého předpokladu M2, řešení soustavy A*X = Bse tedy zefektivní po převodu na ekvivalentní soustavu inv(M1)*A*inv(M2)*Y = inv(M1)*B, kde X = inv(M1)*Y.
35
5 Popis praktické části bakalářské práce Praktická část bakalářské práce je napsaná v programovacím jazyce Java. Je rozdělena na dvě třídy. V třídě Algoritmy jsou implementovány všechny algoritmy, kterými jsou zpracovávána vstupní data matic. Třída Hlavni_okno řeší grafickou část program, naplnění matice zadanými hodnotami a volá příslušné algoritmy. V této se vytvoří instance třídy Panel, do které se vykresluje soustava lineárních rovnic.
Obrázek 1 Struktura praktické části BP
5.1 Popis grafického rozhraní programu Grafické rozhraní je rozděleno na 2 části. V horní části aplikace je text field pro zadání řádu matice, tlačítko zobraz, vykreslí rovnice pro zadání jednotlivých koeficientů aij. Combo boxy, kterými uživatel vybere metodu pro řešení soustavy a kritérium podle kterého se v každé iteraci kontroluje, zda se má pokračovat. Další text field načítá hodnotu kritéria. Tlačítko vypočti, po zadání všech parametrů provádí samotný výpočet. Panel, do kterého se vykresluje soustava lineárních rovnic.
Obrázek 2 Grafické rozhraní BP
36
V dolní části aplikace je editor pane, který se plní výsledky iteračních cyklů a label ve kterém se zobrazuje čas výpočtu.
5.2 Popis použitých algoritmů Metoda matPlusMat si alokuje dvourozměrné pole double, které naplní součtem prvků na odpovídajících si pozicí vstupních matic. Nakonec funkce vrací pole, které je výsledkem součtu matic. Metoda matXmat vrací dvourozměrné pole, které si alokovalo v těle funkce a naplnilo hodnotami odpovídajícími součinu matic. Metoda matXvekt vrací jednorozměrné pole odpovídající součinu matice a vektoru. Metoda vektMinusVekt si alokuje jednorozměrné pole double. Naplní ho rozdíly prvků na odpovídajících si pozicích vektorů, vstupujících do metody v argumentu. Nakonec pole vrátí výsledný vektor. Metoda inverse zapíše do reference, která vstupuje do funkce jako druhý argument inverzní matici k matice, která je prvním argumentem funkce. Metody kriterium1, kriterium2 a kriterium3 vypočítávají hodnotu podle zvoleného kritéria. Kriterium1 vrací reálné číslo double, které vypočte podle předpisu pro absolutní chybu. Kritérium2 vrací číslo podle předpisu relativní chyby a kritérium3 podle předpisu residua. Metoda Jacobiho provádí iterační výpočet na základě vstupních argumentů. Provádí jej na základě následující rovnice.
Po provedení sérii výpočtu se funkcí pro výpočet kritéria přiřadí hodnota tolerance aktuální iterace a porovná se s tolerancí, kterou zvolil uživatel. Po ukončení iteračního výpočtu se vrací matice všech iteračních hodnot neznámých x. Metoda JacobiMaticove provádí iterační výpočet stejně jako metoda Jacobiho, jen za použití rozdílné rovnice . Vrací stejně jako metoda předešlá matici všech iteračních hodnot neznámých x Metoda Gauss_Seidelova provádí iterační výpočet podle následující rovnice.
37
Stejně jako předešlé metody pro výpočet iteračních metod, také tato metoda každé iteraci zjišťuje hodnotu tolerance v jednotlivých iteračních cyklech. Po konci výpočtu vrátí matici iteračních hodnot všech výsledků x. Metoda Gauss_SeidelMaticove analogicky k Jacobiovým se chová stejně jako metoda Gauss_Seidelova, jen s rozdílem použití jiné rovnice pro výpočet.
38
6 Srovnání s matematickým softwarem Matlab Matlab je matematický software vyvíjen společností MathWorks. Název Matlab vznikl zkrácením slov matrix laboratoř, což lze volně přeložit jako maticová laboratoř. Tento překlad plně vystihuje strukturu, která je při výpočtech založena na maticích. „Matlab umožňuje počítání s maticemi, vykreslování 2D a 3D grafů funkcí, implementaci algoritmů, počítačovou simulaci, analýzu a prezentaci dat i vytváření aplikací včetně grafického rozhraní“7
6.1 Srovnání časových nároků při výpočtu Srovnání autorem vytvořenou aplikaci a matematickým softwarem Matlab se provede na deseti pokusech výpočtu stejné soustavy lineární rovnice.
Obrázek 3 Kód pro výpočet Jacobiovy metody v Matlabu
7
http://cs.wikipedia.org/wiki/MATLAB
39
Obrázek 4 Kód pro výpočet G-S metody v Matlabu
Tyto kódy byly zvoleny, aby se časová náročnost měřila na přesně stejných úkonech, které se vykonají ve zdrojovém kódu praktické části bakalářské práce.
Tabulka 3 Srovnaní časové náročnosti Jacobiovy metody
40
Jacobiova metoda 70 60 čas[ms]
50 40 30
Aplikace
20
Matlab
10 0 1
2
3
4
5
6
7
8
9
10
pokus
Graf 1 Srovnání časové náročnosti Jacobiovy metody
Tabulka 4 Srovnaní časové náročnosti Jacobiovy metody maticovým způsobem
41
Jacobiova metoda maticově 70 60
čas[ms]
50 40 30
Aplikace
20
Matlab
10 0 1
2
3
4
5
6
7
8
9
10
pokus
Graf 2 Srovnání časové náročnosti Jacobiovy metody maticovým způsobem
Tabulka 5 Srovnaní časové náročnosti G-S metody
42
Gaussova-Seidelova metoda 35 30 čas[ms]
25 20 15
Aplikace
10
Matlab
5 0 1
2
3
4
5
6
7
8
9
10
pokus
Graf 3 Srovnání časové náročnosti G-S metody
Tabulka 6 Srovnaní časové náročnosti G-S metody maticovým způsobem
43
Gaussova-Seidelova metoda maticevě 35 30 čas[ms]
25 20 15
Aplikace
10
Matlab
5 0 1
2
3
4
5
6
7
8
9
10
pokus
Graf 4 Srovnání časové náročnosti G-S metody maticovým způsobem
6.2 Srovnání přesnosti výsledků Přesnost výsledků se provede na totožné matici jako u srovnávání časových nároků na výpočet.
Tabulka 7 Srovnání přesnosti výsledků Jacobiovy metody
44
Tabulka 8 Srovnání přesnosti výsledků G-S metody
Z předchozích tabulek 7 a 8 je zřejmé, že přesnost výsledků praktické části bakalářské práce zaokrouhlených na 4 desetinná místa, je shodná u Jacobiovy metody po 10 iteracích a G-S metoda po 7 iteracích s výsledky získané z Matlabu.
45
Závěr Teoretická část bakalářské práce má za úkol seznámit čtenáře s problematikou řešení systémů lineárních rovnic, jak přímými metodami, tak i iteračními metodami. Při psaní byla snaha popisovat metody matematickými a přesnými definicemi, ale také je popisovat jednoduchým a snadno pochopitelným způsobem. K podtržení snadnosti pochopení se názorně každá metoda ukázala na typových příkladech s jasnými instrukcemi dalších úkonů. Dále zde srovnáváme výslednou aplikaci, navrženou a realizovanou v praktické části bakalářské práce, se softwarem pro vědecké výpočty Matlabem. Příjemným zjištěním bylo, když časové nároky na dobu výpočtu a všech úkonů, které se dějí při výpočtu v praktické části bakalářské práce, jsou rychlejší než ty samé úkony ve zmiňovaném Matlabu. Při srovnávání přesnosti výsledků došlo ke konstatování, že přesnost jednotlivých kořenů soustavy lineárních rovnic je naprosto stejná jako v Matlabu. V praktické části byly následně implementovány algoritmy, které byly popsány v teoretické části. Pro aplikaci byl zvolen programovací jazyk Java, který je považován za jeden z nejpopulárnějších programovacích jazyků na světě pro svou přenositelnost mezi různými platformami.
46
Literatura 2011. Centrum aplikované matematiky. Západočeská univerzita v Plzni. [Online] 2011. [Citace: 6. 5 2012.] http://www.cam.zcu.cz/~danek/Students/2006_LS/soubory/prednasky/NM_prednaska_03. pdf. 2011. Centrum aplikované matematiky. Západočeská univerzita v Plzni. [Online] 2011. [Citace: 26. 4 2012.] http://www.cam.zcu.cz/~danek/Students/2006_LS/soubory/prednasky/NM_prednaska_03. pdf. Dušek, František. 2000. MATLAB a SIMULINK úvod do používání. Pardubice : Univerzita Pardubice, 2000. str. 146. 80-7194-273-1. 2005. Fakulta aplikovaných věd. Západočeská univerzita v Plzni. [Online] 2005. [Citace: 16. 4 2012.] http://dimatia.fav.zcu.cz/2005/vyuka/zm1/soustavy.pdf. Felcman, Jiří. 2012. 4.4 LU rozklad v obecném případě for nm. Scribd. [Online] 2012. [Citace: 16. 4 2012.] http://www.scribd.com/doc/53207894/21/LU-rozklad-v-obecnemp%C5%99ipad%C4%9B. 2011. Gaussova eliminační metoda - Gaussova eliminace. Aristoteles.Cz - matematika, chemie a fyzika online. [Online] 21. 11 2011. [Citace: 2. 5 2012.] http://www.aristoteles.cz/matematika/linearni_algebra/soustavy/gaussova-eliminacnimetoda.php. Horová, Ivana a Zelenka, Jiří. 2008. Numerické metody. Brno : Masarykova univerzita, 2008. str. 299. 978-80-210-3317-7. Koláček, Jan. 2012. Jan Koláček. Masarykova univerzita. [Online] 2012. [Citace: 15. 4 2012.] http://www.math.muni.cz/~kolacek/vyuka/nummet/prime_metody.pdf. 2006. Matematická sekce MFF UK. Univerzita Karlova - UK. [Online] 3. 10 2006. [Citace: 19. 4 2012.] http://www.karlin.mff.cuni.cz/~tuma/2002/NLinalg2.pdf. 2001. MATLAB - Wikipedie. Wikipedie, otevřená encyklopedie. [Online] 2001. [Citace: 1. 5 2012.] http://cs.wikipedia.org/wiki/MATLAB. Rod, Michal. 2008. Michal Rod. ČVUT v Praze. [Online] 2008. [Citace: 25. 4 2012.] http://rod.aspone.cz/files/LU.pdf. 2001. Soustava lineárních rovnic - Wikipedie. Wikipedie, otevřená encyklopedie. [Online] 2001. [Citace: 29. 4 2012.] http://cs.wikipedia.org/wiki/Soustava_line%C3%A1rn%C3%ADch_rovnic.
47
The MathWorks, Inc. 2012. Mathematics :: Function Reference (MATLAB®). MathWorks - MATLAB and Simulink for Technical Computing. [Online] 2012. [Citace: 6. 5 2012.] http://www.mathworks.com/help/techdoc/ref/f16-5872.html#f16-6565. Vitásek, Emil. 1987. Numerické metody. Praha : SNTL-Nakladatelsví technické literatury, 1987. str. 516. 04-009-87.
48
Příloha A – Dokumentace Dokumentace k Aplikaci je přiložena na CD.
49
Příloha B – Zdrojové kódy Zdrojové kódy Aplikaci jsou přiloženy na CD.
50
