Numerické metody

Uˇcebn´ı texty vysok´ych ˇskol NUMERICK´ E METODY doc. RNDr. Libor ˇ Cerm´ak, CSc. RNDr. Rudolf Hlaviˇcka, CSc. © doc. RNDr. Libor ˇ Cerm´ak, CSc., RNDr. Rudolf Hlaviˇcka, CSc., 2025 Tato publikace podl´eh´a licenci Creative Commons Uved’te autora — Neuˇz´ıvejte d´ılo komerˇcnˇe — Nezpracov´avejte 4.0 Mezin´arodn´ı ISBN 978-80-214-6322-6 Tiˇstˇen´e skriptum Numerick´e metody autor˚u doc. RNDr. Libora ˇ Cerm´aka, CSc. a RNDr. Rudolfa Hlaviˇcky, CSc. vyˇslo v roce 2016. Od roku 2024 je snahou Fakulty strojn´ıho inˇzen´yrstv´ı VUT v Brnˇe postupn´e zpˇr´ıstupˇnov´an´ı uˇcebn´ıch text˚u student˚um v elektronick´e verzi v Digit´aln´ı knihovnˇe VUT. Z PDF verze skript z roku 2016 vzniklo toto e-skriptum. ediˇcn´ı rada FSI VUT, ´unor 2025 Obsah 1´ Uvod do problematiky numerick´ych metod 6 1.1 Chyby v numerick´ych v´ypoˇctech . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2 Reprezentace ˇc´ısel v poˇc´ıtaˇci . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.3 Podm´ınˇenost ´uloh a algoritm˚u . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.4 Cviˇcen´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2ˇ Reˇsen´ı soustav line´arn´ıch rovnic 16 2.1 Pˇr´ım´e metody . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.1.1 Gaussova eliminaˇcn´ı metoda . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.1.2 V´ybˇer hlavn´ıho prvku . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.1.3 Vliv zaokrouhlovac´ıch chyb . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.1.4 Podm´ınˇenost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.2 Iteraˇcn´ı metody . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.3 Cviˇcen´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3 Aproximace funkc´ı 41 3.1 Interpolace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.1.1 Interpolace polynomem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.1.2 Interpolaˇcn´ı splajny . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.1.3 Interpolace funkc´ı v´ıce promˇenn´ych . . . . . . . . . . . . . . . . . . . . . . 53 3.2 Metoda nejmenˇs´ıch ˇctverc˚u . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.3 Cviˇcen´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4 Numerick´y v´ypoˇcet derivace a integr´alu 63 4.1 Numerick´e derivov´an´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.2 Richardsonova extrapolace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.3 Numerick´e integrov´an´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.3.1 Z´akladn´ı formule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.3.2 Sloˇzen´e formule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.3.3 Doplˇnuj´ıc´ı poznatky . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.4 Cviˇcen´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5ˇ Reˇsen´ı neline´arn´ıch rovnic 79 5.1 Urˇcen´ı poˇc´ateˇcn´ı aproximace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.2 Zpˇresˇnuj´ıc´ı metody . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.3 Soustavy neline´arn´ıch rovnic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.4 Cviˇcen´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6 Optimalizace 97 6.1 Jednorozmˇern´a minimalizace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 6.2 Minimalizace funkce v´ıce promˇenn´ych . . . . . . . . . . . . . . . . . . . . . . . . 101 6.3 Cviˇcen´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 Literatura 110 4 Pˇredmluva Tato skripta jsou urˇcena student˚um prvn´ıho roˇcn´ıku Fakulty strojn´ıho inˇzen´yrstv´ı VUT v Brnˇe pro studium pˇredmˇetu Numerick´e metody I. Skripta jsou vˇenov´ana tradiˇcn´ım t´emat˚um numerick´e matematiky: zdroj˚um chyb a jejich ˇs´ıˇren´ı, ˇreˇsen´ı soustav line´arn´ıch rovnic, aproximaci funkc´ı, numerick´emu derivov´an´ı a integrov´an´ı, ˇreˇsen´ı neline´arn´ıch rovnic a minimalizaci funkc´ı. Pˇredpokl´ad´a se, ˇze ˇcten´aˇr skript m´a z´akladn´ı znalosti z line´arn´ı algebry, z diferenci´aln´ıho a integr´aln´ıho poˇctu funkc´ı jedn´e promˇenn´e a z programov´an´ı. Funkce v´ıce promˇenn´ych se v tomto textu vyskytuj´ı jen v odstavc´ıch 3.1.3, 4.3.3, 5.3 a 6.2 a k pochopen´ı zde prob´ıran´e l´atky postaˇc´ı znalosti pr˚ubˇeˇznˇe z´ısk´avan´e v paralelnˇe prob´ıhaj´ıc´ım kurzu MATEMATIKA II, vˇenovan´em pˇrev´aˇznˇe diferenci´aln´ımu a integr´aln´ımu poˇctu funkc´ı v´ıce promˇenn´ych. Do skripta jsme ke kaˇzd´emu t´ematu zaˇradili klasick´e metody, standardnˇe uv´adˇen´e v kaˇzd´em ´uvodn´ım kurzu numerick´e matematiky, kter´e slouˇz´ı pˇredevˇs´ım k pochopen´ı a ilustraci dan´e problematiky. Klasick´e metody jsou dnes jiˇz ˇcasto pˇrekon´any efektivnˇejˇs´ımi postupy. Pot´ıˇz s modern´ımi, v souˇcasnosti pouˇz´ıvan´ymi algoritmy, je vˇsak v tom, ˇze b´yvaj´ı ˇcasto pomˇernˇe komplikovan´e, takˇze porozumˇet jim neb´yv´a snadn´e. Pˇresto se alespoˇn o nˇekter´ych modern´ıch algoritmech v tomto textu struˇcnˇe zmiˇnujeme. Vˇetˇsinu metod jsme se pokusili zd˚uvodnit, precizn´ımu dokazov´an´ı jsme se vˇsak ´umyslnˇe vyh´ybali. U nˇekter´ych tvrzen´ı uv´ad´ıme liter´arn´ı zdroj, v nˇemˇz lze naj´ıt odpov´ıdaj´ıc´ı vysvˇetlen´ı. Jde-li vˇsak o tvrzen´ı natolik bˇeˇzn´e, ˇze ho lze naj´ıt prakticky v jak´ekoliv z´akladn´ı uˇcebnici numerick´ych metod, pak liter´arn´ı zdroj neuv´ad´ıme. ˇ Reˇsen´e pˇr´ıklady jsme volili tak, aby pom´ahaly pochopit a pouˇz´ıvat nˇekdy ponˇekud tˇeˇzko straviteln´e formule a postupy. Skripta obsahuj´ı celou ˇradu algoritm˚u, ˇreˇsen´ych pˇr´ıklad˚u a ´uloh k samostatn´emu procviˇcen´ı. Vˇetˇsina ´uloh se d´a jen tˇeˇzko zvl´adnout bez vyuˇzit´ı poˇc´ıtaˇce. Proto se pˇredpokl´ad´a, ˇze studenti si nˇekter´e algoritmy sami v MATLABu naprogramuj´ı. Pˇri zpracov´an´ı skript jsme vych´azeli z osvˇedˇcen´ych uˇcebnic numerick´e matematiky, jak´ymi jsou napˇr. knihy [3], [22], ale tak´e z modernˇejˇs´ıch knih [17], [7], [11] a [15]. Zvl´aˇstˇe poslednˇe citovan´a kniha Molerova, kter´a je volnˇe ke staˇzen´ı na Internetu, byla pro n´as velkou inspirac´ı. Pokud jde o ˇcesk´e zdroje, nejv´ıce podnˇet˚u jsme ˇcerpali z knih [19], [14], [8]. Modern´ı ˇcesky psan´a monografie numerick´ych metod v souˇcasnosti nen´ı k dispozici. Kromˇe v´yˇse uveden´ych knih lze vˇsak v ˇceˇstinˇe ˇc´ıst tak´e [20], vybran´e kapitoly vˇenovan´e numerick´ym metod´am v [21] a [23]. Douf´ame, ˇze se V´am, mil´y ˇcten´aˇri, numerick´e metody zal´ıb´ı, a ˇze je dok´aˇzete pozdˇeji efektivnˇe vyuˇz´ıt pˇri ˇreˇsen´ı konkr´etn´ıch technick´ych probl´em˚u. Mˇejte ale pros´ım poˇr´ad na zˇreteli, ˇze to, co se o numerick´ych metod´ach v tomto textu dozv´ıte, je opravdu jen naprost´y z´aklad. Zbytek je uˇz na V´as, zdroj˚u informac´ı je dostatek, zejm´ena d´ıky Internetu. Z´akladn´ı orientace v numerick´ych metod´ach V´am pak umoˇzn´ı kvalifikovanˇe pouˇz´ıvat nepˇrebern´e mnoˇzstv´ı program˚u, at ’uˇz volnˇe dostupn´ych na Internetu nebo program˚u komerˇcn´ıch. A to se m˚uˇze hodit. Brno, listopad 2016 Libor ˇ Cerm´ak Rudolf Hlaviˇcka 5 1. ´ Uvod do problematiky numerick´ych metod Pˇri ˇreˇsen´ı probl´em˚u re´aln´eho svˇeta se st´ale ˇcastˇeji setk´av´ame s potˇrebou popsat zkoumanou skuteˇcnost pomoc´ı vˇerohodn´eho matematick´eho modelu a ten pak uspokojivˇe vyˇreˇsit. ˇ Zijeme v dobˇe poˇc´ıtaˇc˚u a tak je pˇrirozen´e, ˇze k realizaci matematick´eho modelu poˇc´ıtaˇc vyuˇzijeme. Poˇc´ıtaˇce um´ı pracovat velmi rychle s informacemi k´odovan´ymi pomoc´ı ˇc´ısel. A pr´avˇe zde je m´ısto pro numerickou matematiku (v angliˇctinˇe numerical analysis ) jakoˇzto vˇedn´ı discipl´ınu, kter´a vyv´ıj´ı a analyzuje metody, jejichˇz technologick´ym j´adrem jsou manipulace s ˇc´ısly. V posledn´ıch letech se v anglicky psan´e literatuˇre m´ısto term´ınu numerical analysis st´ale ˇcastˇeji pouˇz´ıv´a term´ın scientific computing (odpov´ıdaj´ıc´ı ˇcesk´y term´ın n´am bohuˇzel nen´ı zn´am). Kdyˇz chceme metodami numerick´e matematiky vyˇreˇsit dan´y probl´em popsan´y obecn´ym matematick´ym modelem, mus´ıme takov´y model nejdˇr´ıve digitalizovat, to jest formulovat ho ve tvaru numerick´e ´ulohy, jej´ıˇz vstupn´ı i v´ystupn´ı data jsou ˇc´ısla. Numerick´a metoda je postup ˇreˇsen´ı numerick´e ´ulohy. Pˇresn´y popis krok˚u realizuj´ıc´ıch numerickou metodu oznaˇcujeme jako algoritmus numerick´e metody. Lze ho vyj´adˇrit jako posloupnost akc´ı (provediteln´ych na poˇc´ıtaˇci), kter´e k dan´emu (pˇresnˇe specifikovan´emu koneˇcn´emu) souboru vstupn´ıch ˇc´ısel jednoznaˇcnˇe pˇriˇrad´ı odpov´ıdaj´ıc´ı (pˇresnˇe specifikovan´y koneˇcn´y) soubor v´ystupn´ıch ˇc´ısel. Pˇr´ıprava rozs´ahl´ych soubor˚u vstupn´ıch dat b´yv´a oznaˇcov´ana jako preprocessing. Rozsah souboru v´ysledn´ych ´udaj˚u je ˇcasto ohromn´y, pro ˇclovˇeka nestraviteln´y, a proto je tˇreba v´ysledky vhodnˇe zpˇr´ıstupnit tak, aby je zadavatel v´ypoˇctu byl v˚ubec schopen vyhodnotit. Metod´am, kter´e to prov´adˇej´ı, se ˇr´ık´a postprocessing. Jednou z forem postprocessingu je vizualizace v´ysledk˚u. Jako pˇr´ıklad probl´emu ze ˇzivota uvaˇzujme pˇredpovˇed’poˇcas´ı. Pohyb vzduchu v atmosf´eˇre dovedeme alespoˇn pˇribliˇznˇe popsat pomoc´ı soustav parci´aln´ıch diferenci´aln´ıch rovnic a vhodn´ych doplˇnuj´ıc´ıch podm´ınek. Metodami numerick´e matematiky dok´aˇzeme tyto rovnice pˇribliˇznˇe ˇreˇsit. Potˇrebn´a vstupn´ı data se z´ısk´avaj´ı pomoc´ı druˇzic a pozemn´ıch meteorologick´ych stanoviˇst ’. V´ysledky numerick´ych v´ypoˇct˚u zpracovan´e do animovan´ych meteorologick´ych map pak sledujeme v televizn´ı pˇredpovˇedi poˇcas´ı. Pˇri ˇreˇsen´ı re´aln´ych probl´em˚u t´emˇeˇr nikdy nez´ısk´ame pˇresn´e ˇreˇsen´ı, mus´ıme se spokojit jen s ˇreˇsen´ım pˇribliˇzn´ym, kter´e je zat´ıˇzeno chybami. Naˇs´ım c´ılem je organizovat v´ypoˇcet tak, aby celkov´a chyba byla co nejmenˇs´ı. Pˇredevˇs´ım se mus´ıme vyvarovat hrub´ych lidsk´ych chyb, kter´e vypl´yvaj´ı z nepochopen´ı probl´emu a z nepozornosti nebo nedbalosti ˇclovˇeka pˇri jeho ˇreˇsen´ı. Chyba matematick´eho modelu. Pˇri vytv´aˇren´ı matematick´eho modelu re´aln´eho probl´emu prov´ad´ıme vˇzdy jist´e idealizace. Rozd´ıl mezi ˇreˇsen´ım idealizovan´eho probl´emu a ˇreˇsen´ım probl´emu re´aln´eho naz´yv´ame chybou matematick´eho modelu. Do t´eto kategorie chyb zahrnujeme tak´e chyby ve vstupn´ıch ´udaj´ıch. Pˇr´ıklad. M´ame urˇcit povrch zemsk´eho pl´aˇstˇe. K v´ypoˇctu pouˇzijeme vzorec S= 4πr2pro povrch koule o polomˇeru r. Chyba modelu spoˇc´ıv´a v pˇredpokladu, ˇze Zemˇe je koule. Chyba numerick´e metody. Jestliˇze k ˇreˇsen´ı (numerick´e) ´ulohy pouˇzijeme numerickou metodu, kter´a n´am neposkytne pˇresn´e (teoretick´e) ˇreˇsen´ı dan´e ´ulohy, pak chybu, kter´e se dopust´ıme, naz´yv´ame chybou numerick´e metody. D˚uleˇzitou souˇc´ast´ı n´avrhu numerick´e 6 metody je odhad chyby numerick´e metody. Pˇr´ıklad. M´ame spoˇc´ıtat hodnotu funkce sin 1 seˇcten´ım koneˇcn´eho poˇctu ˇclen˚u Taylorovy ˇrady sin x=x−x3 3! +x5 5! −x7 7! +x9 9! −···+ (−1)nx2n+1 (2n+ 1)! +··· pro x= 1. Je zn´amo, ˇze seˇcten´ım prvn´ıch tˇr´ı ˇclen˚u ˇrady se dopust´ıme chyby velikosti nejv´yˇse 1/7!, obecnˇe seˇcten´ım prvn´ıch nˇclen˚u se dopust´ıme chyby nejv´yˇse 1/(2n+ 1)!. Zaokrouhlovac´ı chyby. Pˇri pr´aci na poˇc´ıtaˇci m˚uˇzeme k reprezentaci ˇc´ısel pouˇz´ıt jen koneˇcn´y poˇcet cifer. Pracujeme proto s pˇribliˇzn´ymi hodnotami ˇc´ısel, kter´e dostaneme zaokrouhlen´ım pˇresn´ych hodnot. Zaokrouhlovac´ı chyby vznikaj´ı uˇz pˇri vkl´ad´an´ı dat do poˇc´ıtaˇce, dalˇs´ı pak vznikaj´ı pˇri ˇc´ıseln´ych v´ypoˇctech. Pˇri ˇspatnˇe organizovan´em v´ypoˇctu m˚uˇze doj´ıt v d˚usledku nahromadˇen´ı zaokrouhlovac´ıch chyb k naprost´emu znehodnocen´ı v´ysledku, viz pˇr´ıklad 1.7. Pˇr´ıklad. ˇ C´ıslo πneum´ıme do poˇc´ıtaˇce vloˇzit pˇresnˇe. Tak´e v´ysledek operace, pˇri n´ıˇz ˇc´ıslo 2 dˇel´ıme ˇc´ıslem 3, nezobraz´ıme na standardn´ım poˇc´ıtaˇci pracuj´ıc´ım s bin´arn´ımi ˇc´ısly pˇresnˇe. Je tˇreba m´ıt na pamˇeti, ˇze pˇri ˇreˇsen´ı re´aln´eho probl´emu vystupuj´ı obvykle vˇsechny chyby souˇcasnˇe. 1.1. Chyby v numerick´ych v´ypoˇctech Absolutn´ı a relativn´ı chyba. Ve v´ypoˇctech jsme ˇcasto nuceni nahradit pˇresn´e ˇc´ıslo x pˇribliˇzn´ym ˇc´ıslem ˜x.ˇ C´ıslo ˜xpotom naz´yv´ame aproximac´ı ˇc´ısla x. Rozd´ıl ˜x−x= ∆x naz´yv´ame absolutn´ı chybou aproximace ˜xa ˇc´ıslo ∆x x=˜x−x x, x 6= 0 , naz´yv´ame relativn´ı chybou aproximace ˜x. Pro |∆x| ≤ εse pouˇz´ıv´a tak´e symbolick´y z´apis ˜x=x±εa m´ın´ı se t´ım, ˇze x−ε≤˜x≤x+ε. Podobnˇe se pro |∆x/x| ≤ δpouˇz´ıv´a z´apis ˜x=x(1 ±δ). Absolutn´ı hodnota relativn´ı chyby se ˇcasto uv´ad´ı v procentech. Nyn´ı posoud´ıme chybu, kter´e se dopust´ıme pˇri v´ypoˇctu hodnoty f(x1, x2,...,xn) funkce f, kdyˇz pˇresn´e hodnoty xinahrad´ıme pˇribliˇzn´ymi hodnotami ˜xi=xi+ ∆xi. Z Taylorova rozvoje f(˜ x) okolo bodu xdostaneme f(˜ x) = f(x) + n X i=1 ∆xi ∂f(x) ∂xi +1 2 n X i,j=1 ∆xi∆xj ∂2f(x) ∂xi∂xj +···. Povaˇzujeme-li souˇciny chyb ∆xi∆xjza mal´e, m´ame pro absolutn´ı chybu |∆f(x)|:= |f(˜ x)−f(x)|. = n X i=1 ∂f(x) ∂xi ∆xi≤ n X i=1  ∂f(x) ∂xi·|∆xi|(1.1) 7 a pro chybu relativn´ı  ∆f(x) f(x). = n X i=1 xi f(x) ∂f(x) ∂xi ∆xi xi≤ n X i=1  xi f(x) ∂f(x) ∂xi· ∆xi xi.(1.2) Pˇri praktick´ych odhadech se hodnota funkce fa hodnoty jej´ıch derivac´ı ∂f/∂xina prav´ych stran´ach pˇribliˇzn´ych nerovnost´ı (1.2) a (1.1) poˇc´ıtaj´ı v bodˇe ˜ x. Chyby z´akladn´ıch aritmetick´ych operac´ı. Zvol´ıme-li f(x, y) = x±y, dostaneme pro absolutn´ı a relativn´ı chybu souˇctu a rozd´ılu ∆(x±y). = ∆x±∆y , ∆(x±y) x±y . =x x±y ∆x x±y x±y ∆y y.(1.3) Pro vyj´adˇren´ı chyby souˇcinu vol´ıme f(x, y) = xy a obdrˇz´ıme ∆(xy). =y∆x+x∆y , ∆(xy) xy . =∆x x+∆y y(1.4) a pro chybu pod´ılu dostaneme volbou f(x, y) = x/y ∆x y. =1 y∆x−x y2∆y , ∆(x/y) x/y . =∆x x−∆y y.(1.5) Vˇsimnˇete si, ˇze relativn´ı chyba souˇctu resp. rozd´ılu m˚uˇze b´yt v´yraznˇe vˇetˇs´ı neˇz relativn´ı chyby operand˚u v pˇr´ıpadˇe, kdyˇz |x±y|je podstatnˇe menˇs´ı neˇz |x|nebo |y|. Pˇri dˇelen´ı mal´ym ˇc´ıslem je (d´ıky druh´e mocninˇe yve jmenovateli) v´yznamn´a chyba absolutn´ı. Platn´e dekadick´e cifry. Necht ’˜xje aproximace ˇc´ısla x, kterou zapiˇsme v mocninn´em dekadick´em rozvoji jako ˜x=±[d1·10e+d2·10e−1+···+dk·10e+1−k+dk+1 ·10e−k+...], d16= 0 . ˇ Rekneme, ˇze k−t´a dekadick´a cifra dkaproximace ˜xje platn´a, jestliˇze |˜x−x| ≤ 5·10e−k,(1.6) tj. kdyˇz se ˜xliˇs´ı od xnejv´yˇse o 5 jednotek ˇr´adu pˇr´ısluˇsn´eho n´asleduj´ıc´ı cifˇre. Plat´ı-li nerovnost (1.6) pro k≤p, ale pro k=p+ 1 uˇz neplat´ı, ˇr´ık´ame, ˇze ˜xm´a pplatn´ych cifer. ˇ C´ıslo ˜x=±d1,d2d3. . . dp·10e, kter´e m´a vˇsech pcifer platn´ych, je spr´avnˇe zaokrouhlenou hodnotou ˇc´ısla x. Platn´a desetinn´a m´ısta. ˇ Rekneme, ˇze aproximace ˜xˇc´ısla xm´a k-t´e desetinn´e m´ısto platn´e, jestliˇze |˜x−x| ≤ 5·10−k−1,(1.7) tj. kdyˇz se ˜xliˇs´ı od xnejv´yˇse o 5 jednotek ˇr´adu pˇr´ısluˇsn´eho n´asleduj´ıc´ımu desetinn´emu m´ıstu. Plat´ı-li nerovnost (1.7) pro k≤p, ale pro k=p+ 1 uˇz neplat´ı, ˇr´ık´ame, ˇze ˜x m´a pplatn´ych desetinn´ych m´ıst. Ve spr´avnˇe zaokrouhlen´em ˇc´ısle je tedy kaˇzd´e desetinn´e m´ısto platn´e. V n´asleduj´ıc´ı tabulce uv´ad´ıme nˇekolik pˇr´ıklad˚u: 8 x˜xplatn´e cifry platn´a desetinn´a m´ısta 284 290 1 − −45,8472 −45,798 3 1 100,002 99,9973 4 2 99,9973 100,002 5 2 −0,003728 −0,0041 1 3 1,841 ·10−62,5·10−60 5 Pˇri odeˇc´ıt´an´ı dvou bl´ızk´ych ˇc´ısel doch´az´ı ke ztr´atˇe platn´ych cifer, jak o tom svˇedˇc´ı Pˇr´ıklad 1.1. Je-li x= 4,998949 ·101,˜x= 4,999 ·101,|∆x|= 5,10 ·10−4, ∆x x. = 1,020 ·10−5, y= 5,001848 ·101,˜y= 5,002 ·101,|∆y|= 1,52 ·10−3, ∆y y. = 3,039 ·10−5, pak pro rozd´ıly z=y−x, ˜z= ˜y−˜xdost´av´ame z= 2,899 ·10−2,˜z= 3 ·10−2,|∆z|= 1,01 ·10−3, ∆z z. = 3,484 ·10−2, takˇze ˜zm´a jen jednu platnou cifru, zat´ımco ˜xi ˜ymaj´ı ˇctyˇri platn´e cifry.  Pˇr´ıklad 1.2. Necht ’x= 1,3262 ±5·10−5,y=−6,5347 ±5·10−5,z= 13,235 ±5·10−4. M´ame urˇcit aproximaci funkˇcn´ı hodnoty f=xy/z, absolutn´ı a relativn´ı chybu a poˇcet platn´ych cifer v´ysledku. Spoˇcteme ˜ f= ˜x˜y/˜z=−6,548031 ...·10−1.Podle (1.1) pak pˇribliˇznˇe plat´ı  ∆f ˜ f≤ ˜y ˜z∆x+ ˜x ˜z∆y+ ˜x˜y ˜z2∆z ˜x˜y ˜z −1 = ∆x ˜x+ ∆y ˜y+ ∆z ˜z. = 8,31 ·10−5. Odtud |∆f|. = 8,31 ·10−5·|˜ f|. = 5,44 ·10−5<5·10−1−3,takˇze (se tˇremi platn´ymi ciframi) f=−0,6548 ±0,0001.  1.2. Reprezentace ˇc´ısel v poˇc´ıtaˇci Re´aln´a ˇc´ısla jsou v poˇc´ıtaˇc´ıch reprezentov´ana v syst´emu ˇc´ısel s pohyblivou ˇr´adovou ˇc´arkou (v angliˇctinˇe floating point numbers). Z´akladn´ı myˇslenka je podobn´a semilogaritmick´emu z´apisu (v angliˇctinˇe scientific notation), v nˇemˇz napˇr. ˇc´ıslo 245700 p´ıˇseme jako 2,457·105a ˇc´ıslo 0,0005768 jako 5,768·10−4. V tomto form´atu se desetinn´a ˇc´arka pohybuje (v doslovn´em pˇrekladu plave) v z´avislosti na dekadick´em exponentu. Form´alnˇe lze syst´em Fnormalizovan´ych ˇc´ısel pohybliv´e ˇr´adov´e ˇc´arky charakterizovat ˇctyˇrmi cel´ymi ˇc´ısly: βz´aklad ˇc´ıseln´e soustavy (β≥2), ppˇresnost (p≥1), [L , U] rozsah exponentu (L < 0< U). 9 2. ˇ Reˇsen´ı soustav line´arn´ıch rovnic Jednou z nejˇcastˇeji se vyskytuj´ıc´ıch ´uloh v´ypoˇcetn´ı praxe je ´uloha vyˇreˇsit soustavu line´arn´ıch rovnic. Takov´e soustavy b´yvaj´ı ˇcasto velmi rozs´ahl´e, souˇcasn´a v´ypoˇcetn´ı technika umoˇzˇnuje v pˇrijateln´ych ˇcasech vyˇreˇsit soustavy s nˇekolika mili´ony nezn´am´ych. Metody ˇreˇsen´ı dˇel´ıme na pˇr´ım´e a iteraˇcn´ı. Pˇr´ım´e metody jsou takov´e metody, kter´e dodaj´ı v koneˇcn´em poˇctu krok˚u pˇresn´e ˇreˇsen´ı za pˇredpokladu, ˇze v´ypoˇcet prob´ıh´a bez zaokrouhlovac´ıch chyb, tedy zcela pˇresnˇe. Iteraˇcn´ı metody poskytnou jen ˇreˇsen´ı pˇribliˇzn´e. To ale v˚ubec nevad´ı, pokud je pˇribliˇzn´e ˇreˇsen´ı dostateˇcnˇe dobrou aproximac´ı ˇreˇsen´ı pˇresn´eho. Poˇcet krok˚u iteraˇcn´ı metody z´avis´ı na poˇzadovan´e pˇresnosti. Budeme se tedy zab´yvat ˇreˇsen´ım soustavy line´arn´ıch rovnic a11 x1+a12 x2+··· +a1nxn=b1, a21 x1+a22 x2+··· +a2nxn=b2, . . .. . . an1x1+an2x2+··· +ann xn=bn. (2.1) Soustavu (2.1) m˚uˇzeme ps´at ve tvaru n X j=1 aijxj=bi, i = 1,2,...,n, (2.2) nebo v maticov´em tvaru Ax =b,(2.3) kde A=     a11 a12 ··· a1n a21 a22 ··· a2n . . .. . .. . .. . . an1an2··· ann     ,x=     x1 x2 . . . xn     ,b=     b1 b2 . . . bn     . Matici Anaz´yv´ame matic´ı soustavy,bje vektor prav´e strany axvektor nezn´am´ych. Budeme pˇredpokl´adat, ˇze matice soustavy je regul´arn´ı, takˇze ˇreˇsen´a soustava m´a jedin´e ˇreˇsen´ı. 2.1. Pˇr´ım´e metody 2.1.1. Gaussova eliminaˇcn´ı metoda Z´akladn´ı pˇr´ımou metodou ˇreˇsen´ı soustav line´arn´ıch rovnic je Gaussova eliminaˇcn´ı metoda, struˇcnˇe GEM. Skl´ad´a se ze dvou ˇc´ast´ı. V pˇr´ım´em chodu GEM se soustava (2.1) pˇrevede na ekvivalentn´ı soustavu Ux =c,(2.4) 16 kde Uje tzv. horn´ı troj´uheln´ıkov´a matice, coˇz je matice, kter´a m´a pod hlavn´ı diagon´alou vˇsechny prvky nulov´e, tj. U={uij}n i,j=1 auij = 0 pro i > j, U=          u11 u12 u13 ··· u1,n−1u1n 0u22 u23 ··· u2,n−1u2n 0 0 u33 ··· u3,n−1u3n . . .. . .. . .. . .. . .. . . 0 0 ··· 0un−1,n−1un−1,n 0 0 0 ··· 0unn          . Ve zpˇetn´em chodu se pak ˇreˇs´ı soustava (2.4). Protoˇze Aje regul´arn´ı, je tak´e Uregul´arn´ı, coˇz znamen´a, ˇze diagon´aln´ı prvky uii 6= 0, i= 1,2,...,n. D´ıky tomu vypoˇcteme z posledn´ı rovnice xn, z pˇredposledn´ı xn−1atd. aˇz nakonec z prvn´ı rovnice vypoˇcteme x1. Pˇr´ım´y chod GEM. Pro usnadnˇen´ı popisu pˇr´ım´eho chodu GEM poloˇz´ıme A(0) =A, b(0) =b, prvky matice matice A(0) oznaˇc´ıme a(0) ij ≡aij a prvky vektoru b(0) oznaˇc´ıme b(0) i≡bi. Pˇr´ım´y chod GEM popisuje n´asleduj´ıc´ı algoritmus GEMz (z´akladn´ı, bez v´ybˇeru hlavn´ıho prvku): for k:= 1 to n−1do begin A(k):= A(k−1) ;b(k):= b(k−1) ; for i:= k+ 1 to ndo begin mik := a(k) ik /a(k) kk ; for j:= k+ 1 to ndo a(k) ij := a(k) ij −mika(k) kj ; b(k) i:= b(k) i−mikb(k) k; end end Pˇr´ım´y chod m´a n−1 krok˚u. V k-t´em kroku se soustava rovnic A(k−1)x=b(k−1) transformuje na soustavu A(k)x=b(k). Prvn´ıch krovnic se uˇz nemˇen´ı. Tato skuteˇcnost je v algoritmu GEMz vyj´adˇrena pˇr´ıkazy A(k):= A(k−1) ab(k):= b(k−1). Smyslem transformace je vylouˇcit nezn´amou xkz rovnic i > k, tj. vynulovat poddiagon´aln´ı koeficienty v k-t´em sloupci matice A(k). Dos´ahneme toho tak, ˇze od i-t´e rovnice odeˇcteme mik n´asobek k-t´e rovnice. Multiplik´atory mik musej´ı zajistit, aby v pozici (i, k) matice A(k)vznikla nula: a(k) ik −mika(k) kk = 0 =⇒mik =a(k) ik /a(k) kk . ˇ C´ıslo a(k) kk se naz´yv´a hlavn´ı prvek nebo tak´e pivot. Pˇri v´ypoˇctu multiplik´atoru mik m˚uˇze algoritmus GEMz zhavarovat v pˇr´ıpadˇe, ˇze a(k) kk = 0. Tomuto probl´emu bychom se mohli vyhnout, kdybychom k-tou rovnici prohodili s nˇekterou z dalˇs´ıch rovnic, kter´a m´a u promˇenn´e xknenulov´y koeficient. Postup zaloˇzen´y na t´eto myˇslence se naz´yv´a GEM s v´ybˇerem hlavn´ıho prvku. Podrobnˇe se j´ım budeme zab´yvat v n´asleduj´ıc´ım odstavci. GEMz je tedy algoritmus GEM bez v´ybˇeru hlavn´ıho prvku. 17 V tomto odstavci budeme pˇredpokl´adat, ˇze Aje takov´a matice soustavy, pro kterou jsou vˇsechny hlavn´ı prvky a(k) kk nenulov´e. Programov´an´ı. Prvky matic A(k)uchov´av´ame v dvourozmˇern´em poli Aa prvky vektor˚u b(k)v jednorozmˇern´em poli b. Pˇr´ıkazy A(k):= A(k−1) ab(k):= b(k−1) se proto ve skuteˇcnosti neprov´adˇej´ı. Protoˇze v pozici (i, k) pole Avznikne nula, lze prvek A(i,k) vyuˇz´ıt pro uskladnˇen´ı multiplik´atoru mik. LU rozklad. Po ukonˇcen´ı pˇr´ım´eho chodu je horn´ı troj´uheln´ıkov´a matice Uv rovnici (2.4) urˇcena diagon´aln´ımi a naddiagon´aln´ımi prvky matice A(n−1), tj. uij := 0 pro j= 1,2,...,i−1, a(n−1) ij pro j=i, i + 1,...,n, i= 1,2,...,n. (2.5) Vektor cv (2.4) je transformovanou pravou stranou b(n−1), tj. ci:= b(n−1) i, i = 1,2,...,n. (2.6) Multiplik´atory mij z pˇr´ım´eho chodu um´ıst´ıme do doln´ı troj´uheln´ıkov´e matice L={ℓij}n i,j=1, ℓij = 0 pro j > i, L=          ℓ11 0 0 ··· 0 0 ℓ21 ℓ22 0··· 0 0 ℓ31 ℓ32 ℓ33 ··· 0 0 . . .. . .. . .. . .. . .. . . ℓn−1,1ℓn−1,2ℓn−1,3. . . ℓn−1,n−10 ℓn1ℓn2ℓn3··· ℓn,n−1ℓnn          , definovan´e pˇredpisem ℓik :=    0 pro i= 1,2,...,k−1, 1 pro i=k , mik pro i=k+ 1, k + 2,...,n, k= 1,2,...,n. (2.7) D´a se uk´azat, ˇze plat´ı A=LU .(2.8) Vyj´adˇren´ı matice Ajako souˇcinu doln´ı troj´uheln´ıkov´e matice La horn´ı troj´uheln´ıkov´e matice Use naz´yv´a LU rozklad matice A. Ten je moˇzn´e pouˇz´ıt k pozdˇejˇs´ımu ˇreˇsen´ı soustavy rovnic se stejnou matic´ı soustavy A, avˇsak s jinou pravou stranou. To je uˇziteˇcn´e zejm´ena pˇri ˇreˇsen´ı posloupnosti ´uloh Axi=bi, kdy se nov´a prav´a strana bim˚uˇze sestavit aˇz pot´e, co se vyˇreˇsily pˇredchoz´ı soustavy Axk=bkpro k < i. Ukaˇzme si, jak lze soustavu LUx =befektivnˇe vyˇreˇsit. Kdyˇz si oznaˇc´ıme Ux =y, vid´ıme, ˇze yje ˇreˇsen´ı soustavy Ly =b. Urˇc´ıme tedy nejdˇr´ıve yjako ˇreˇsen´ı soustavy Ly =ba pak xjako ˇreˇsen´ı soustavy Ux =y, tj. Ly =b,Ux =y.(2.9) 18 Zˇrejmˇe y=b(n−1) je transformovan´a prav´a strana z´ıskan´a algoritmem GEMz. Soustavu Ly =bvyˇreˇs´ıme snadno, z prvn´ı rovnice vypoˇc´ıt´ame y1, ze druh´e rovnice y2atd. aˇz nakonec z posledn´ı rovnice vypoˇc´ıt´ame yn. Soustavu Ux =yˇreˇs´ıme pozp´atku, tj. z posledn´ı rovnice vypoˇcteme xn, z pˇredposledn´ı xn−1atd. aˇz nakonec z prvn´ı rovnice vypoˇcteme x1. Pˇri ˇreˇsen´ı soustav rovnic b´yv´a LU rozklad oznaˇcov´an tak´e jako eliminace nebo pˇr´ım´y chod a v´ypoˇcet ˇreˇsen´ı podle (2.9) b´yv´a oznaˇcov´an jako zpˇetn´y chod. Kdy lze algoritmus GEMz pouˇz´ıt? Jak jsme jiˇz uvedli, slab´ym m´ıstem algoritmu GEMz m˚uˇze b´yt v´ypoˇcet multiplik´atoru mik, nebot ’obecnˇe nelze vylouˇcit, ˇze v pr˚ubˇehu eliminace vznikne a(k) kk = 0. V aplikac´ıch se vˇsak pomˇernˇe ˇcasto ˇreˇs´ı soustavy rovnic, pro kter´e nulov´y pivot v algoritmu GEMz vzniknout nem˚uˇze. Abychom takov´e soustavy mohli popsat, zavedeme si nˇekolik nov´ych pojm˚u. ˇ Rekneme, ˇze matice A={aij}n i,j=1 je ryze diagon´alnˇe dominantn´ı, jestliˇze |aii|> n X j= 1 j6=i |aij|, i = 1,2,...,n, (2.10) nebo-li slovy, v kaˇzd´em ˇr´adku je absolutn´ı hodnota diagon´aln´ıho prvku vˇetˇs´ı neˇz souˇcet absolutn´ıch hodnot zb´yvaj´ıc´ıch prvk˚u tohoto ˇr´adku. I kdyˇz matice Asoustavy rovnic Ax =bdiagon´alnˇe dominantn´ı nen´ı, lze nˇekdy vhodn´ym pˇreskl´ad´an´ım rovnic doc´ılit toho, ˇze matice ˆ Atakto vznikl´e ekvivalentn´ı soustavy rovnic ˆ Ax =ˆ buˇz diagon´alnˇe dominantn´ı je. V aplikac´ıch se tak´e pomˇernˇe ˇcasto setk´av´ame s tzv. pozitivnˇe definitn´ımi maticemi. Takov´e matice lze specifikovat pomoc´ı ˇrady navz´ajem ekvivalentn´ıch definic. Jednu z nich si ted’uvedeme: ˇrekneme, ˇze matice A={aij}n i,j=1 je pozitivnˇe definitn´ı, jestliˇze je symetrick´a a pro kaˇzd´y nenulov´y sloupcov´y vektor x= (x1, x2,...,xn)Tplat´ı xTAx = n X i,j=1 xiaijxj>0.(2.11) Ovˇeˇrit pˇr´ımo tuto podm´ınku nen´ı snadn´e. Je-li vˇsak Aregul´arn´ı, pak z (2.11) okamˇzitˇe plyne, ˇze ATAje pozitivnˇe definitn´ı. (Dokaˇzte to!) Vyn´asob´ıme-li tedy soustavu rovnic Ax =bzleva matic´ı AT, dostaneme ekvivalentn´ı soustavu ATAx =ATbs pozitivnˇe definitn´ı matic´ı soustavy. Tento postup se vˇsak pro praktick´e ˇreˇsen´ı soustav rovnic nehod´ı (operace ATAvyˇzaduje velk´y objem v´ypoˇct˚u, u iteraˇcn´ıch metod se nav´ıc v´yznamnˇe zhorˇsuje rychlost konvergence). Pˇri ˇreˇsen´ı konkr´etn´ıch praktick´ych ´uloh b´yv´a obvykle uˇz pˇredem zn´amo (z povahy ˇreˇsen´eho probl´emu a ze zp˚usobu jeho diskretizace), zda matice vznikaj´ıc´ıch soustav line´arn´ıch rovnic jsou (resp. nejsou) pozitivnˇe definitn´ı. Uved’me si vˇsak pˇresto alespoˇn jednu ˇcasto uv´adˇenou (nutnou a postaˇcuj´ıc´ı) podm´ınku pozitivn´ı definitnosti, zn´amou jako Sylvesterovo krit´erium. ˇ Ctvercov´a symetrick´a matice A={aij}n i,j=1 je pozitivnˇe definitn´ı, pr´avˇe kdyˇz jsou kladn´e determinanty vˇsech hlavn´ıch rohov´ych submatic {aij}k i,j=1, 19 k= 1,2,...,n, tj. kdyˇz plat´ı a11 >0,a11 a12 a21 a22>0, a11 a12 a13 a21 a22 a23 a31 a32 a33>0, . . . ,  a11 . . . a1n . . .. . . an1. . . ann  >0. D´a se uk´azat, ˇze algoritmus GEMz lze pouˇz´ıt pro ˇreˇsen´ı soustav, jejichˇz matice je bud’to ryze diagon´alnˇe dominantn´ı nebo pozitivnˇe definitn´ı.´ Uspˇeˇsn´e pouˇzit´ı algoritmu GEMz lze zaruˇcit tak´e pro dalˇs´ı typy matic, kter´e se pˇri ˇreˇsen´ı praktick´ych ´uloh ˇcasto vyskytuj´ı (viz napˇr. [4], [13]). V´ypoˇctov´a n´aroˇcnost GEM. Pˇr´ım´y chod GEM vyˇzaduje 1 3n3+O(n2) operac´ı n´asobic´ıch (tj. n´asoben´ı nebo dˇelen´ı) a 1 3n3+O(n2) operac´ı sˇc´ıtac´ıch (tj. sˇc´ıt´an´ı nebo odeˇc´ıt´an´ı). Symbolem O(n2) jsme pˇritom vyj´adˇrili ˇr´adovˇe m´enˇe v´yznamn´y poˇcet operac´ı ˇr´adu n2 (tvaru α2n2+α1n+α0, kde α2,α1,α0jsou ˇc´ısla nez´avisl´a na n). ˇ Clen 1 3n3souvis´ı s transformac´ı matice soustavy. Poˇcet operac´ı souvisej´ıc´ıch s transformac´ı prav´e strany je o ˇr´ad niˇzˇs´ı a je tedy zahrnut do ˇclenu O(n2). Zpˇetn´y chod GEM je v´ypoˇcetnˇe podstatnˇe m´enˇe n´aroˇcn´y. ˇ Reˇsen´ı soustavy rovnic s troj´uheln´ıkovou matic´ı vyˇzaduje 1 2n2+O(n) operac´ı n´asobic´ıch a 1 2n2+O(n) operac´ı sˇc´ıtac´ıch. Pˇritom O(n) reprezentuje poˇcet operac´ı ˇr´adu n(tvaru α1n+α0, kde α1,α0jsou ˇc´ısla nez´avisl´a na n). GEM zpˇetn´y chod, tj. v´ypoˇcet xze soustavy (2.4), proto vyˇzaduje 1 2n2+O(n) operac´ı a LU zpˇetn´y chod, tj. v´ypoˇcet yaxze soustav (2.9), vyˇzaduje dvojn´asobn´y poˇcet operac´ı, tj. n2+O(n). Pro velk´y poˇcet rovnic, tj. pro velk´e n, proto m˚uˇzeme tvrdit, ˇze eliminace vyˇzaduje pˇribliˇznˇe 1 3n3operac´ı a GEM resp. LU zpˇetn´y chod pˇribliˇznˇe 1 2n2resp. n2operac´ı (n´asobic´ıch a stejnˇe tak sˇc´ıtac´ıch). Cholesk´eho rozklad. Pozitivnˇe definitn´ı matici Alze vyj´adˇrit ve tvaru A=LLT,(2.12) kde Lje doln´ı troj´uheln´ıkov´a matice, jej´ıˇz nenulov´e prvky jsou postupnˇe pro k= 1,2,...,n urˇceny pˇredpisem ℓkk =v u u takk − k−1 X j=1 ℓ2 kj ,(2.13) ℓik =1 ℓkk aik − k−1 X j=1 ℓijℓkj!, i =k+ 1, k + 2,...,n. Soustavu rovnic ˇreˇs´ıme podle (2.9) pro U=LT. Vyj´adˇren´ı matice Ave tvaru (2.12) se naz´yv´a Cholesk´eho rozklad matice A. Cholesk´eho rozklad vyˇzaduje pˇribliˇznˇe poloviˇcn´ı v´ypoˇctov´e n´aklady oproti obecn´emu LU rozkladu, tedy pˇribliˇznˇe 1 6n3operac´ı n´asobic´ıch a zhruba stejn´y poˇcet operac´ı sˇc´ıtac´ıch (v´ypoˇcet odmocnin nem´a na celkov´y poˇcet operac´ı podstatn´y vliv). Cholesk´eho algoritmus (2.13) lze pouˇz´ıt k efektivn´ımu posouzen´ı pozitivn´ı definitnosti matice A: je-li Asymetrick´a a plat´ı-li akk −Pk−1 j=1 ℓ2 kj >0 pro k= 1,2,...,n, pak Aje pozitivnˇe definitn´ı. V MATLABu lze pro Cholesk´eho rozklad pouˇz´ıt funkci chol. 20 2.1.2. V´ybˇer hlavn´ıho prvku Zaˇcneme pˇr´ıkladem. Pˇr´ıklad 2.1. M´ame vyˇreˇsit soustavu rovnic   10 −7 0 −3 2,099 6 5−1,1 4,8   x1 x2 x3 =  7 3,901 5,9  na hypotetick´em poˇc´ıtaˇci, kter´y pracuje v dekadick´e soustavˇe s pˇetim´ıstnou mantisou. Pˇresn´e ˇreˇsen´ı je x=  0 −1 1 . V prvn´ım kroku eliminujeme poddiagon´aln´ı prvky v prvn´ım sloupci a dostaneme   10 −7 0 0−0,001 6 0 2,4 4,8   x1 x2 x3 =  7 6,001 2,4 . Prvek v pozici (2,2) je ve srovn´an´ı s ostatn´ımi prvky matice mal´y. Pˇresto pokraˇcujme v eliminaci. V dalˇs´ım kroku je tˇreba ke tˇret´ımu ˇr´adku pˇriˇc´ıst ˇr´adek druh´y n´asoben´y 2400: (4,8 + 6 ·2400) ·x3= 2,4 + 6,001 ·2400 . Na lev´e stranˇe je koeficient 4,8+6·2400 = 14404,8 zaokrouhlen na 14405. Na prav´e stranˇe v´ysledek n´asoben´ı 6,001 ·2400 = 14402,4 nelze zobrazit pˇresnˇe, mus´ı b´yt zaokrouhlen na 14402. K tomu se pak pˇriˇcte 2,4 a znovu dojde k zaokrouhlen´ı. Posledn´ı rovnice tak nabude tvaru 14405 x3= 14404 . Zpˇetn´y chod zaˇcne v´ypoˇctem x3=14404 14 405 . = 0,99993 . Pˇresn´y v´ysledek je x3= 1. Zd´a se, ˇze chyba nen´ı nijak v´aˇzn´a. Bohuˇzel, x2je tˇreba urˇcit z rovnice −0,001 x2+ 6 ·0,99993 = 6,001, coˇz d´av´a, po zaokrouhlen´ı 6 ·0,99993 . = 5,9996, x2=0,0014 −0,001 =−1,4. Nakonec vypoˇcteme x1z prvn´ı rovnice 10x1−7·(−1,4) = 7 21 a dostaneme x1=−0,28. M´ısto pˇresn´eho ˇreˇsen´ı xjsme dostali pˇribliˇzn´e ˇreˇsen´ı ˜ x= −0,28 −1,4 0,99993 . Kde vznikl probl´em? Nedoˇslo k ˇz´adn´emu hromadˇen´ı chyb zp˚usoben´emu prov´adˇen´ım tis´ıc˚u operac´ı. Matice soustavy nen´ı bl´ızk´a matici singul´arn´ı. Pot´ıˇz je jinde, p˚usob´ı ji mal´y pivot ve druh´em kroku eliminace. T´ım vznikne multiplik´ator −2400 a v d˚usledku toho m´a posledn´ı rovnice koeficienty zhruba 1000 kr´at vˇetˇs´ı neˇz koeficienty p˚uvodn´ı rovnice. Zaokrouhlovac´ı chyby, kter´e jsou mal´e vzhledem k tˇemto velk´ym koeficient˚um, jsou nepˇrijateln´e pro koeficienty p˚uvodn´ı matice a tak´e pro samotn´e ˇreˇsen´ı. Snadno se provˇeˇr´ı, ˇze kdyˇz druhou a tˇret´ı rovnici prohod´ıme, nevzniknou ˇz´adn´e velk´e multiplik´atory a v´ysledek je zcela pˇresn´y. Ukazuje se, ˇze to plat´ı obecnˇe: jestliˇze jsou absolutn´ı hodnoty multiplik´ator˚u menˇs´ı nebo nejv´yˇse rovny 1, pak je numericky spoˇcten´e ˇreˇsen´ı vyhovuj´ıc´ı.  ˇ C´asteˇcn´y v´ybˇer hlavn´ıho prvku je modifikace GEM zajiˇst ’uj´ıc´ı, aby absolutn´ı hodnota multiplik´ator˚u byla menˇs´ı nebo rovna jedn´e. V k-t´em kroku eliminace se jako pivot vyb´ır´a prvek s nejvˇetˇs´ı absolutn´ı hodnotou v zat´ım neeliminovan´e ˇc´asti k-t´eho sloupce matice A(k−1), tj. mezi prvky a(k−1) ik pro i≥k. Necht ’tedy rje takov´y ˇr´adkov´y index, pro kter´y |a(k−1) rk |= max k≤i≤n|a(k−1) ik |.(2.14) Pak prohod´ıme k-tou a r-tou rovnici. Ze soustavy rovnic A(k−1)x=b(k−1) tak dostaneme soustavu A(k)x=b(k), pˇriˇcemˇz A(k)z´ısk´ame prohozen´ım k-t´eho a r-t´eho ˇr´adku matice A(k−1) a podobnˇe b(k)z´ısk´ame prohozen´ım k-t´eho a r-t´eho prvku vektoru b(k−1). Poddiagon´aln´ı prvky v k-t´em sloupci matice A(k)eliminujeme stejnˇe jako v algoritmu GEMz. ˇ Reˇsen´ı soustavy line´arn´ıch rovnic s ˇc´asteˇcn´ym v´ybˇerem hlavn´ıch prvk˚u v MATLABu dostaneme pomoc´ı pˇr´ıkazu x=A\b.                  a(k−1) 11 a(k−1) 12 ··· a(k−1) 1k ··· a(k−1) 1n 0a(k−1) 22 ··· a(k−1) 2k ··· a(k−1) 2n . . . . . . . . . . . . . . . . . . 0 0 ··· a(k−1) kk ··· a(k−1) kn . . . . . . . . . . . . . . . . . . 0 0 ··· a(k−1) rk ··· a(k−1) rn . . . . . . . . . . . . . . . . . . 0 0 ··· a(k−1) nk ··· a(k−1) nn                                   x1 x2 . . . xk . . . xr . . . xn                  =                  b(k−1) 1 b(k−1) 2. . . b(k−1) k . . . b(k−1) r . . . b(k−1) n                  Obr. 2.1: GEM s ˇc´asteˇcn´ym v´ybˇerem hlavn´ıho prvku (v krouˇzku) ´ Upln´y v´ybˇer hlavn´ıho prvku je postup, kter´y m˚uˇze absolutn´ı hodnoty multiplik´ator˚u zmenˇsit jeˇstˇe v´yraznˇeji. Dociluje se toho t´ım, ˇze v k-t´em kroku eliminace se jako pivot vyb´ır´a prvek s nejvˇetˇs´ı absolutn´ı hodnotou v dosud neeliminovan´e ˇc´asti matice A(k−1), tj. 22 v ˇr´adc´ıch i≥ka sloupc´ıch j≥k. Necht ’tedy rje ˇr´adkov´y a ssloupcov´y index vybran´y tak, ˇze |a(k−1) rs |= max k≤i,j≤n|a(k−1) ij |.(2.15) Pak prohod´ıme k-tou a r-tou rovnici a k-tou a s-tou nezn´amou. Ze soustavy rovnic                        a(k−1) 11 a(k−1) 12 ··· a(k−1) 1k ··· a(k−1) 1r···a(k−1) 1s···a(k−1) 1n 0a(k−1) 22 ··· a(k−1) 2k ··· a(k−1) 2r···a(k−1) 2s···a(k−1) 2n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 ··· a(k−1) kk ··· a(k−1) kr ···a(k−1) ks ···a(k−1) kn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 ··· a(k−1) rk ··· a(k−1) rr ···a(k−1) rs ···a(k−1) rn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 ··· a(k−1) sk ··· a(k−1) sr ···a(k−1) ss ···a(k−1) sn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 ··· a(k−1) nk ··· a(k−1) nr ···a(k−1) ns ···a(k−1) nn                                               x(k−1) 1 x(k−1) 2 . . . x(k−1) k . . . x(k−1) r . . . x(k−1) s . . . x(k−1) n                        =                        b(k−1) 1 b(k−1) 2. . . b(k−1) k . . . b(k−1) r . . . b(k−1) s . . . b(k−1) n                        Obr. 2.2: GEM s ´upln´ym v´ybˇerem hlavn´ıho prvku (v krouˇzku) A(k−1)x(k−1) =b(k−1) dostaneme soustavu A(k)x(k)=b(k), pˇriˇcemˇz A(k)z´ısk´ame prohozen´ım k-t´eho a r-t´eho ˇr´adku a k-t´eho a s-t´eho sloupce matice A(k−1),b(k)z´ısk´ame prohozen´ım k-t´eho a r-t´eho prvku vektoru b(k−1) ax(k)z´ısk´ame prohozen´ım k-t´eho as-t´eho prvku vektoru x(k−1). Pˇritom x(0) =xje p˚uvodn´ı vektor nezn´am´ych. Prohazov´an´ı promˇenn´ych registrujeme ve vektoru p= (p1, p2,...,pn)T. Na poˇc´atku poloˇz´ıme pi=ia pˇri kaˇzd´em prohozen´ı promˇenn´ych prohod´ıme tak´e odpov´ıdaj´ıc´ı prvky vektoru p. Po ukonˇcen´ı pˇr´ım´eho chodu je i-t´a sloˇzka vektoru x(n−1) rovna pi-t´e sloˇzce p˚uvodn´ıho vektoru x. Ve zpˇetn´em chodu vypoˇcteme x(n−1) a pomoc´ı vektoru purˇc´ıme x. ˇ C´asteˇcn´y nebo ´upln´y v´ybˇer hlavn´ıho prvku? Vˇetˇsinou se pouˇz´ıv´a jen ˇc´asteˇcn´y v´ybˇer hlavn´ıho prvku. Praxe i teorie potvrzuje, ˇze i ˇc´asteˇcn´y v´ybˇer hlavn´ıch prvk˚u staˇc´ı k tomu, aby zaokrouhlovac´ı chyby z˚ustaly dostateˇcnˇe mal´e a neznehodnotily v´ysledn´e ˇreˇsen´ı. Dalˇs´ım d˚uvodem, proˇc ˇc´asteˇcn´emu v´ybˇeru hlavn´ıch prvk˚u d´av´ame pˇrednost, je to, ˇze jeho realizace vyˇzaduje v´yraznˇe menˇs´ı poˇcet operac´ı neˇz v´ybˇer ´upln´y. LU rozklad s ˇc´asteˇcn´ym v´ybˇerem hlavn´ıho prvku je standardn´ı rutina dostupn´a v kaˇzd´e knihovnˇe program˚u pro numerick´e ˇreˇsen´ı ´uloh line´arn´ı algebry (v MATLABu viz funkce lu). Vstupn´ım parametrem je matice A. V´ystupn´ı parametry jsou tˇri: doln´ı troj- ´uheln´ıkov´a matice L(s jedniˇckami na hlavn´ı diagon´ale), horn´ı troj´uheln´ıkov´a matice U a tzv. permutaˇcn´ı matice P, pˇriˇcemˇz LU =PA .(2.16) Pˇritom permutaˇcn´ı matice je takov´a matice, kter´a vznikne z jednotkov´e matice nˇejak´ym proh´azen´ım jej´ıch ˇr´adk˚u. N´asleduje popis algoritmu pro LU rozklad matice As ˇc´asteˇcn´ym v´ybˇerem hlavn´ıch prvk˚u. 23 Algoritmus LUp (s ˇc´asteˇcn´ym v´ybˇerem hlavn´ıho prvku) 1) Poloˇz´ıme A(0) =A,P(0) =I. 2) Postupnˇe pro k= 1,2,...,n−1 prov´ad´ıme 2a) Urˇc´ıme index rpivotn´ıho ˇr´adku, viz (2.14). 2b) Prohod´ıme ˇr´adky karv matici A(k−1) a takto z´ıskanou matici oznaˇc´ıme jako A(k). Prohod´ıme rovnˇeˇz ˇr´adky karv matici P(k−1) a takto z´ıskanou matici oznaˇc´ıme jako P(k). 2c) Uprav´ıme matici A(k)tak, ˇze eliminujeme poddiagon´aln´ı prvky v k-t´em sloupci, tj. postupnˇe pro i=k+ 1, k + 2,...,n poˇc´ıt´ame mik := a(k) ik /a(k) kk , a(k) ij := a(k) ij −mika(k) kj pro j=k+ 1, k + 2,...,n, a(k) ik := mik . Do anulovan´ych pozic (i, k) matice A(k)tedy ukl´ad´ame multiplik´atory mik. 3) Doln´ı troj´uheln´ıkovou matici Lsestroj´ıme tak, ˇze do hlavn´ı diagon´aly d´ame jedniˇcky a poddiagon´aln´ı prvky pˇrevezmeme z v´ysledn´e matice A(n−1). Horn´ı troj´uheln´ıkovou matici Udostaneme z diagon´aln´ıch a naddiagon´aln´ıch prvk˚u v´ysledn´e matice A(n−1). Permutaˇcn´ı matice Pje rovna v´ysledn´e permutaˇcn´ı matici P(n−1). Po z´ısk´an´ı matic L,UaPvyˇreˇs´ıme uˇz soustavu rovnic Ax =bsnadno. Zˇrejmˇe PAx =Pb =⇒LUx =Pb . Proto nejdˇr´ıve urˇc´ıme vektor z=Pb, tj. prvky vektoru bprav´e strany proh´az´ıme stejnˇe, jako jsme prohazovali ˇr´adky matic A(k−1) aP(k−1) v algoritmu LUp. Oznaˇc´ıme-li Ux =y, vid´ıme, ˇze yje ˇreˇsen´ı soustavy Ly =z. Vyˇreˇsen´ım t´eto soustavy dostaneme y. Zb´yv´a jeˇstˇe vyˇreˇsit soustavu Ux =ya ˇreˇsen´ı xje nalezeno. Shrneme-li to, poˇc´ıt´ame postupnˇe z,yaxz rovnic z=Pb ,Ly =z,Ux =y.(2.17) Uchov´avat historii prohazov´an´ı ˇr´adk˚u v matici Pje zˇrejm´e pl´ytv´an´ı pamˇet ’ov´ym m´ıstem, vˇzdyt ’z celkov´eho poˇctu n2prvk˚u matice Pje jich pouze nnenulov´ych. Proto m´ısto s matic´ı Pstaˇc´ı pracovat jen s vektorem ˇr´adkov´ych permutac´ı p. Na zaˇc´atku poloˇz´ıme p(0) = (1,2,...,n)Ta ve zbytku algoritmu pak prohazujeme prvky vektoru p(k−1). Matice Pbyla zapotˇreb´ı jen pro sestaven´ı vektoru z. To ale dok´aˇzeme s pomoc´ı vektoru ptak´e: do i-t´e sloˇzky vektoru zvloˇz´ıme pi-tou sloˇzku vektoru b, postupnˇe pro i= 1,2,...,n. Pˇr´ıklad 2.2. Provedeme LU rozklad matice A=    −0,4−0,95 −0,4−7,34 0,5−0,3 2,15 −2,45 −2 4 1 −3 −1 5,5 2,5 3,5    . 24 V prvn´ım sloupci najdeme jako pivota ˇc´ıslo −2 ve tˇret´ım ˇr´adku. Proto prohod´ıme prvn´ı a tˇret´ı ˇr´adek. Pak eliminujeme poddiagon´aln´ı prvky v prvn´ım sloupci a na jejich m´ısta zap´ıˇseme pouˇzit´e multiplik´atory. Prohozen´ı ˇr´adk˚u vyznaˇc´ıme tak´e v permutaˇcn´ım vektoru. Tak dostaneme A(1) =    −2 4 1 −3 −0,25 0,7 2,4−3,2 0,2−1,75 −0,6−6,74 0,5 3,5 2 5    ,p(1) =    3 2 1 4     Prvn´ı ˇr´adek se uˇz mˇenit nebude. Ve druh´em kroku najdeme pivota ve druh´em sloupci. Je to ˇc´ıslo 3,5 ve ˇctvrt´em ˇr´adku. Proto prohod´ıme druh´y a ˇctvrt´y ˇr´adek jak v matici A(1) tak ve vektoru p(1). Pak eliminujeme prvky v pozic´ıch (3,2) a (4,2) a na jejich m´ısta zap´ıˇseme pouˇzit´e multiplik´atory. V´ysledkem je A(2) =    −2 4 1 −3 0,5 3,5 2 5 0,2−0,5 0,4−4,24 −0,25 0,2 2 −4,2    ,p(2) =    3 4 1 2    . Prvn´ı dva ˇr´adky uˇz z˚ustanou bez zmˇeny. Pivot ve tˇret´ım sloupci je ˇc´ıslo 2 ve ˇctvrt´em ˇr´adku. Prohod´ıme tedy tˇret´ı a ˇctvrt´y ˇr´adek v matici A(2) i ve vektoru p(2). Pak eliminujeme prvek v pozici (4,3) a na jeho m´ısto vloˇz´ıme pouˇzit´y multiplik´ator. Tak dostaneme A(3) =    −2 4 1 −3 0,5 3,5 2 5 −0,25 0,2 2 −4,2 0,2−0,5 0,2−3,4    ,p(3) =    3 4 2 1    . Dostali jsme tedy L=    1 0 0 0 0,5 1 0 0 −0,25 0,2 1 0 0,2−0,5 0,2 1    ,U=    −2 4 1 −3 0 3,5 2 5 0 0 2 −4,2 0 0 0 −3,4    . Permutaˇcn´ı vektor p=p(3). Pokud bychom chtˇeli vytvoˇrit permutaˇcn´ı matici P, staˇc´ı vz´ıt jednotkovou matici a pˇreuspoˇr´adat ji tak, ˇze p˚uvodnˇe pi-t´y ˇr´adek se stane ˇr´adkem i-t´ym. Kdyˇz to provedeme, dostaneme pro permutaˇcn´ı vektor p=    3 4 2 1    permutaˇcn´ı matici P=    0010 0001 0100 1000    . Snadno se ovˇeˇr´ı, ˇze LU =PA. Ukaˇzme si jeˇstˇe ˇreˇsen´ı soustavy rovnic pro volenou pravou stranu. Zvolme tˇreba b=    −13,14 2,15 9 27,5    ,pak z=    9 27,5 2,15 −13,14    , 25 Symbolem Ojsme si pˇritom oznaˇcili nulovou matici, tj. matici, jej´ıˇz vˇsechny prvky jsou rovny nule. Poznamenejme, ˇze obecnˇe je norma ˇctvercov´e matice definov´ana jako re´aln´a funkce splˇnuj´ıc´ı jen prvn´ı ˇctyˇri z v´yˇse uveden´ych podm´ınek. Pˇr´ıklad 2.4. Soustava rovnic Ax =b, kde A=1 10 10 101,b=11 111,m´a ˇreˇsen´ı x=1 1. Pouˇzijeme l∞-normu a spoˇcteme kbk∞= 111, kxk∞= 1 . Kdyˇz pravou stranu zmˇen´ıme na ˜ b=11,11 110,89,dostaneme ˇreˇsen´ı ˜ x=13,21 −0,21. Oznaˇc´ıme-li ∆b=˜ b−b, ∆x=˜ x−x, pak k∆bk∞= 0,11 a k∆xk∞= 12,21. Vid´ıme, ˇze pomˇernˇe mal´a zmˇena prav´e strany zcela zmˇenila ˇreˇsen´ı. Relativn´ı zmˇeny jsou k∆bk∞ kbk∞ = 9,909 ·10−4,k∆xk∞ kxk∞ = 12,21 . Podle (2.21) m˚uˇzeme odhadnout ˇc´ıslo podm´ınˇenosti κ(A)≥12,21 9,909 ·10−4= 12321. Ve skuteˇcnosti je ba ∆bzvoleno tak, ˇze κ(A) = 12321. To se snadno ovˇeˇr´ı, nebot ’ A−1=101 −10 −10 1,takˇze kA−1k∞= 111 = kAk∞aκ(A) = 1112= 12321. Ukaˇzme si jeˇstˇe, ˇze vztah (2.21) plat´ı jako rovnost: k˜ x−xk∞ kxk∞ =12,21 1= 11120,11 111 =κ(A)k∆bk∞ kbk∞ = 111 ·k∆bk∞= 111 ·krk∞, kde r=b−A˜ xje reziduum.  K urˇcen´ı κ(A) potˇrebujeme zn´at kA−1k. Avˇsak v´ypoˇcet A−1vyˇzaduje pˇribliˇznˇe tˇrikr´at tolik pr´ace jako cel´e ˇreˇsen´ı soustavy rovnic. Naˇstˇest´ı pˇresnou hodnotu κ(A) obvykle nepotˇrebujeme, vystaˇc´ıme s dostateˇcnˇe dobr´ym odhadem κ(A). Spolehliv´e a pomˇernˇe velmi rychl´e odhady ˇc´ısla podm´ınˇenosti matic patˇr´ı v souˇcasn´e dobˇe ke standardn´ımu vybaven´ı program˚u pro ˇreˇsen´ı soustav line´arn´ıch rovnic. Jestliˇze program zjist´ı, ˇze ˇc´ıslo podm´ınˇenosti je pˇr´ıliˇs velk´e, vyd´a varov´an´ı nebo dokonce v´ypoˇcet pˇreruˇs´ı. Shrnut´ı. Soustava line´arn´ıch rovnic je dobˇre (ˇspatnˇe) podm´ınˇen´a, pr´avˇe kdyˇz je matice soustavy dobˇre (ˇspatnˇe) podm´ınˇen´a. Podm´ınˇenost matice soustavy Amˇeˇr´ıme pomoc´ı ˇc´ısla podm´ınˇenosti κ(A)≥1. Je-li ˇc´ıslo κ(A) mal´e, je matice Adobˇre podm´ınˇen´a. V opaˇcn´em pˇr´ıpadˇe, tj. kdyˇz κ(A)≫1, je matice Aˇspatnˇe podm´ınˇen´a. ˇ Spatnˇe podm´ınˇenou soustavu rovnic lze obvykle jen velmi obt´ıˇznˇe ˇreˇsit. Pomoci m˚uˇze v´ypoˇcet s v´ıcem´ıstnou mantisou 32 (je vhodn´e pouˇz´ıt dvojn´asobnou nebo jeˇstˇe vˇetˇs´ı pˇresnost). Existuj´ı vˇsak v´yjimky: je-li napˇr´ıklad Adiagon´aln´ı matice, ve kter´e aii = 10i, pak je κ(A) = 10n−1, coˇz je pro velk´e nvelk´e ˇc´ıslo, a pˇresto ˇreˇsen´ı xi= 10−ibiz´ısk´ame bez probl´em˚u pro libovolnˇe velk´y poˇcet rovnic. Pˇredpokl´adejme, ˇze matice soustavy je dobˇre podm´ınˇen´a. Pak je GEM s ˇc´asteˇcn´ym (nebo ´upln´ym) v´ybˇerem hlavn´ıho prvku dobˇre podm´ınˇen´y algoritmus: protoˇze velikost multiplik´ator˚u nepˇresahuje jedniˇcku, vznikaj´ıc´ı zaokrouhlovac´ı chyby se dalˇs´ım v´ypoˇctem nezesiluj´ı. Kdyˇz naopak v´ybˇer hlavn´ıch prvk˚u neprov´ad´ıme, m˚uˇzeme dostat multiplik´atory, jejichˇz absolutn´ı hodnota je vˇetˇs´ı neˇz jedna, coˇz m´a za n´asledek zvˇetˇsov´an´ı dˇr´ıve vznikl´ych zaokrouhlovac´ıch chyb. GEM bez v´ybˇeru hlavn´ıch prvk˚u je tedy obecnˇe ˇspatnˇe podm´ınˇen´y algoritmus. V´yjimku z tohoto pravidla pˇredstavuje ˇreˇsen´ı soustav se speci´aln´ı matic´ı soustavy, napˇr. kdyˇz je matice soustavy ostˇre diagon´alnˇe dominantn´ı nebo pozitivnˇe definitn´ı, pak je i GEM bez v´ybˇeru hlavn´ıch prvk˚u dobˇre podm´ınˇen´y algoritmus. 2.2. Iteraˇcn´ı metody Mnoho praktick´ych probl´em˚u vyˇzaduje ˇreˇsen´ı rozs´ahl´ych soustav line´arn´ıch rovnic Ax =b, v nichˇz matice Aje naˇstˇest´ı ˇr´ıdk´a, tj. m´a relativnˇe m´alo nenulov´ych prvk˚u. Standardn´ı eliminaˇcn´ı metody studovan´e v pˇredchoz´ı kapitole 2.1 nejsou pro ˇreˇsen´ı takov´ych soustav vhodn´e, nebot ’v pr˚ubˇehu eliminace doch´az´ı postupnˇe k zaplˇnov´an´ı p˚uvodnˇe nenulov´ych pozic v matici soustavy, coˇz vede k velk´ym n´arok˚um na poˇcet aritmetick´ych operac´ı a klade tak´e vysok´e n´aroky na pamˇet ’poˇc´ıtaˇce. To je d˚uvod, proˇc se pro ˇreˇsen´ı takov´ych soustav pouˇz´ıvaj´ı iteraˇcn´ı metody. Zvol´ı se poˇc´ateˇcn´ı vektor x0a generuje se posloupnost vektor˚u x0→x1→x2..., kter´a konverguje k hledan´emu ˇreˇsen´ı x. Spoleˇcn´ym rysem vˇsech iteraˇcn´ıch metod je fakt, ˇze kaˇzd´y jednotliv´y iteraˇcn´ı krok xk→xk+1 vyˇzaduje objem v´ypoˇct˚u srovnateln´y s n´asoben´ım matice A vektorem, coˇz je pro ˇr´ıdk´e matice objem nevelk´y (pokud je v kaˇzd´em ˇr´adku matice A ˇr´adu nnejv´yˇse mnenulov´ych prvk˚u, jde o nm operac´ı n´asoben´ı a sˇc´ıt´an´ı). Pˇrijateln´y objem v´ypoˇct˚u lze proto dos´ahnout i pro pomˇernˇe velk´y poˇcet iterac´ı. Na obhajobu pˇr´ım´ych metod je vˇsak tˇreba dodat, ˇze pro soustavy s ˇr´ıdk´ymi maticemi existuj´ı tak´e velmi efektivn´ı algoritmy eliminaˇcn´ıho typu. Pˇresto, pro extr´emnˇe rozs´ahl´e soustavy rovnic se speci´aln´ı strukturou matice soustavy jsou vhodnˇe zvolen´e iteraˇcn´ı metody efektivnˇejˇs´ı a jsou ˇcasto jedinou prakticky realizovatelnou metodou ˇreˇsen´ı. Konvergence. Vˇetˇsina klasick´ych iteraˇcn´ıch metod vych´az´ı z rozkladu matice soustavy A=M−N,kde Mje regul´arn´ı matice. Pak je posloupnost xkdefinov´ana pˇredpisem Mxk+1 =Nxk+b,(2.26) pˇriˇcemˇz poˇc´ateˇcn´ı aproximace x0je dan´a. ˇ Rekneme, ˇze iteraˇcn´ı metoda konverguje, a p´ıˇseme xk→x, kdyˇz ˇc´ıseln´a posloupnost kxk−xk → 0. Oznaˇcme ek=xk−xchybu v k-t´e iteraci. Protoˇze Mx =Nx +b, dostaneme M(xk+1 −x) = N(xk−x) nebo-li ek+1 =M−1Nek. 33 Oznaˇc´ıme-li T=M−1N, pak pomoc´ı (2.24) dost´av´ame kek+1k ≤ kTk·kekk ≤ kTk2·kek−1k ≤ ··· ≤ kTkk+1 ·ke0k. Konvergence iteraˇcn´ı metody z libovoln´eho startovac´ıho vektoru proto jistˇe nastane, kdyˇz kTk<1,kde T=M−1Nje iteraˇcn´ı matice. (2.27) Podm´ınka (2.27) se neovˇeˇruje snadno, a proto si u konkr´etn´ıch metod uvedeme jin´e postaˇcuj´ıc´ı podm´ınky konvergence. Krit´eria pro ukonˇcen´ı iterac´ı. Jde o to, jak rozhodnout, zda xk+1 je uˇz dostateˇcnˇe dobr´a aproximace ˇreˇsen´ı x.ˇ Reˇsen´ı xnezn´ame, takˇze se bez nˇej mus´ıme obej´ıt. Nab´ız´ı se zkoumat velikost zmˇeny xk+1 −xknebo velikost rezidua rk+1 =b−Axk+1. Postupuje se tak, ˇze uˇzivatel zad´a mal´e kladn´e ˇc´ıslo εjako poˇzadovanou pˇresnost a v kaˇzd´em kroku metody se testuje, zda je uˇz splnˇena napˇr. jedna z n´asleduj´ıc´ıch podm´ınek 1. kxk+1 −xkk ≤ εkxkk, 2. krk+1k ≤ ε(kAk·kxk+1k+kbk) , 3. krk+1k ≤ εkr0k. Je-li podm´ınka na ukonˇcen´ı iterac´ı splnˇena, v´ypoˇcet pˇreruˇs´ıme a xk+1 povaˇzujeme za pˇribliˇznou hodnotu ˇreˇsen´ı x. Jacobiova metoda. Pˇredpokl´adejme, ˇze A=L+D+U, kde Dje diagon´aln´ı matice, kter´a m´a stejnou diagon´alu jako A, a kde Lresp. Uje ryze doln´ı resp. horn´ı troj´uheln´ıkov´a ˇc´ast A, tj. D=     a11 0··· 0 0a22 0 . . ..... . . 0 0 ··· ann     , L=     0··· ··· 0 a21 0 0 . . ........ . . an1··· an,n−10     ,U=     0a12 ··· a1n . . ........ . . 0 0 an−1,n 0··· ··· 0     . Nejjednoduˇsˇs´ı rozklad Adostaneme pro M=DaN=−(L+U). Metoda (2.27) je pak tvaru Dxk+1 =b−(L+U)xk(2.28) a je zn´ama jako Jacobiova metoda. Soustava (2.28) s diagon´aln´ı matic´ı se ˇreˇs´ı snadno. Zap´ıˇseme-li (2.28) po sloˇzk´ach (sloˇzky vektoru xkjsou znaˇceny x(k) i, podobnˇe pro xk+1), dostaneme x(k+1) i=1 aii bi− n X j= 1 j6=i aij x(k) j!, i = 1,2,...,n. 34 Anal´yzou vlastnost´ı iteraˇcn´ı matice T=−D−1(L+U) lze dok´azat, ˇze Jacobiova metoda konverguje, kdyˇz Aje ryze diagon´alnˇe dominantn´ı. Gaussova-Seidelova metoda. Vˇsimnˇete si, ˇze Jacobiova metoda pouˇz´ıv´a xkk v´ypoˇctu vˇsech sloˇzek xk+1. Protoˇze (alespoˇn na s´eriov´ych poˇc´ıtaˇc´ıch) prvky vektoru xk+1 poˇc´ıt´ame postupnˇe jeden za druh´ym, vznikl pˇrirozen´y n´apad vyuˇz´ıt ihned ty sloˇzky xk+1, kter´e jsou uˇz k dispozici. Tak dost´av´ame Gaussovu-Seidelovu metodu: x(k+1) i=1 aii bi− i−1 X j=1 aij x(k+1) j− n X j=i+1 aij x(k) j!, i = 1,2,...,n. Vyj´adˇr´ıme-li tuto metodu v maticov´em tvaru, m´ame (D+L)xk+1 =b−Uxk. Je dok´az´ano, ˇze Gaussova-Seidelova metoda konverguje, kdyˇz Aje ryze diagon´alnˇe dominantn´ı nebo pozitivnˇe definitn´ı. Pozn´amky. N´asleduje nˇekolik poznatk˚u o vz´ajemn´em vztahu Jacobiovy a GaussovySeidelovy metody. 1. Konvergence Gaussovy-Seidelovy metody je pro mnoh´e matice Arychlejˇs´ı neˇz konvergence Jacobiovy metody. Tak je tomu tˇreba v pˇr´ıpadˇe, kdyˇz Aje ryze diagon´alnˇe dominantn´ı. 2. Existuj´ı matice, pro kter´e Gaussova-Seidelova metoda konverguje a Jacobiova metoda nekonverguje a naopak, pro kter´e konverguje Jacobiova metoda a GaussovaSeidelova metoda nekonverguje. 3. Jacobiova metoda umoˇzˇnuje paraleln´ı v´ypoˇcet (vˇsechny sloˇzky x(k+1) imohou b´yt poˇc´ıt´any souˇcasnˇe, kaˇzd´a na jin´em procesoru), zat´ımco Gaussova-Seidelova metoda je ze sv´e podstaty sekvenˇcn´ı (x(k+1) ilze vypoˇc´ıtat aˇz po t´e, co byly spoˇcteny vˇsechny sloˇzky x(k+1) jpro j < i). Pro speci´aln´ı typy matic Ajsou vˇsak vypracov´any postupy umoˇzˇnuj´ıc´ı paralelizovat i Gaussovu-Seidelovu metodu. Relaxaˇcn´ı metody. Bezprostˇrednˇe pot´e, co jsme z´akladn´ı metodou (Jacobiovou nebo Gaussovou-Seidelovou) spoˇcetli i-tou sloˇzku x(k+1) i, provedeme jej´ı modifikaci x(k+1) i:= (1 −ω)x(k) i+ωx(k+1) i, kde ω > 0 je tzv. relaxaˇcn´ı parametr. Vol´ıme ho tak, abychom vylepˇsili konvergenci z´akladn´ı metody. Pro ω= 1 dost´av´ame p˚uvodn´ı metodu. Zvol´ıme-li ω < 1, hovoˇr´ıme odoln´ı relaxaci, v pˇr´ıpadˇe ω > 1 jde o horn´ı relaxaci. Efektivn´ı volba relaxaˇcn´ıho parametru ωz´avis´ı na zvolen´e z´akladn´ı metodˇe a na matici soustavy A. Praktick´e zkuˇsenosti potvrzuj´ı, ˇze doln´ı relaxace m˚uˇze zajistit konvergenci v pˇr´ıpadˇe, kdyˇz z´akladn´ı metoda nekonverguje. Vhodnou volbou relaxaˇcn´ıho parametru lze rychlost 35 konvergence p˚uvodn´ı metody podstatnˇe zrychlit. Pro zvolenou metodu a speci´aln´ı tvar matice Ajsou zn´amy vzorce pro optim´aln´ı hodnotu ωopt relaxaˇcn´ıho parametru. Tyto vzorce vˇsak maj´ı v´yznam sp´ıˇse teoretick´y, nebot ’v´ypoˇcet podle nich je pˇr´ıliˇs n´aroˇcn´y. Proto se pracuje s promˇenn´ym relaxaˇcn´ım faktorem, v k-t´e iteraci s ωk, a jeho hodnota se v kaˇzd´e iteraci zpˇresˇnuje tak, aby se postupnˇe bl´ıˇzila k optim´aln´ımu ωopt. Konkr´etn´ı metody lze naj´ıt ve specializovan´e literatuˇre. Relaxace Jacobiovy metody. D´a se uk´azat, ˇze kdyˇz konverguje Jacobiova metoda, tak konverguje tak´e relaxovan´a Jacobiova metoda pro 0 < ω ≤1. Relaxace Gaussovy-Seidelovy metody je v literatuˇre zn´ama jako SOR metoda (podle anglick´eho Successive Over Relaxation). O konvergenci SOR metody m´ame zejm´ena n´asleduj´ıc´ı poznatky: 1. Pokud SOR metoda konverguje, pak je 0 < ω < 2. 2. SOR metoda konverguje, kdyˇz Aje ryze diagon´alnˇe dominantn´ı a 0 < ω ≤1. 3. SOR metoda konverguje, kdyˇz Aje pozitivnˇe definitn´ı a 0 < ω < 2. SOR metoda zapsan´a po sloˇzk´ach je tvaru x(k+1) i=ω aii bi− i−1 X j=1 aij x(k+1) j− n X j=i+1 aij x(k) j!+ (1 −ω)x(k) i, i = 1,2,...,n, a SOR maticovˇe je (ω−1D+L)xk+1 =b−((1 −ω−1)D+U)xk. Kdyˇz ve sloˇzkov´em z´apisu SOR metody na prav´e stranˇe m´ısto x(k+1) jp´ıˇseme x(k) j, dostaneme relaxaci Jacobiovy metody oznaˇcovanou jako JOR metoda, maticovˇe ω−1Dxk+1 =b−(L+ (1 −ω−1)D+U)xk. Pˇr´ıklad 2.5. Velk´e ˇr´ıdk´e matice vznikaj´ı pˇri numerick´em ˇreˇsen´ı parci´aln´ıch diferenci´aln´ıch rovnic. MATLAB m´a ve sv´e galerii pˇr´ıklady takov´ych matic. Pˇr´ıkazem K = gallery(’poisson’,n) vygenerujeme matici, jej´ıˇz strukturu nyn´ı pop´ıˇseme. 1 1Parci´aln´ı diferenci´aln´ı rovnice jsou nezbytn´e pro modelov´an´ı technick´ych probl´em˚u. Pˇri ˇreˇsen´ı Dirichletovy ´ulohy pro Poissonovu rovnici na ˇctverci Ω = (0, l)×(0, l), −∂2u(x, y) ∂x2−∂2u(x, y) ∂y2=f(x, y) v Ω a u(x, y) = g(x, y) na hranici ∂Ω , (funkce fagjsou dan´e, nezn´am´a je funkce u) metodou s´ıt´ı vznik´a soustava line´arn´ıch rovnic s matic´ı soustavy Kuvaˇzovanou v tomto pˇr´ıkladu 2.5. 36 Matice Kje blokovˇe tˇr´ıdiagon´aln´ı, pro devˇet rovnic vypad´a takto               4−1 0 −1 0 0 0 0 0 −1 4 −1 0 −1 0 0 0 0 0−1 4 0 0 −1 0 0 0 −1 0 0 4 −1 0 −1 0 0 0−1 0 −1 4 −1 0 −1 0 0 0 −1 0 −1 4 0 0 −1 000−1 0 0 4 −1 0 0 0 0 0 −1 0 −1 4 −1 00000−1 0 −1 4               . Obecnˇe je Ksloˇzena z ˇctvercov´ych submatic ˇr´adu n(tzv. blok˚u) B,−I,O, kter´e jsou ve ˇctvercov´e matici ˇr´adu n2rozm´ıstˇeny ve tˇrech diagon´al´ach K=           B−I O ··· O O −I B −I··· O O O−I B ··· O O . . .. . .. . ..... . .. . . O O O ··· B−I O O O ··· −I B           . Matice Ije jednotkov´a matice, Oje ˇctvercov´a matice s nulov´ymi prvky a Bje tˇr´ıdiagon´aln´ı matice B=          4−1 0 ··· 0 0 −1 4 −1··· 0 0 0−1 0 ··· 0 0 . . .. . .. . ..... . .. . . 0 0 0 ··· 4−1 0 0 0 ··· −1 4          . Na matici Kse ˇcasto testuje ´uˇcinnost numerick´ych metod pro ˇreˇsen´ı soustav line´arn´ıch rovnic s ˇr´ıdkou matic´ı. Na obr´azku 2.3 je vidˇet z´avislost poˇctu iterac´ı (potˇrebn´ych pro dosaˇzen´ı zvolen´e pˇresnosti) na poˇctu rovnic n2. Pˇri metodˇe SOR bylo zvoleno optim´aln´ı ω. Pozn´amka. Iteraˇcn´ı metody, se kter´ymi jsme se sezn´amili, b´yvaj´ı oznaˇcov´any jako klasick´e iteraˇcn´ı metody. V souˇcasnosti se pouˇz´ıvaj´ı uˇz jen zˇr´ıdka. Do uˇcebn´ıho textu jsme je zaˇradili hlavnˇe proto, ˇze jsou pomˇernˇe jednoduch´e a pˇritom na nich lze dobˇre uk´azat, jak iteraˇcn´ı metody funguj´ı. Na druh´e stranˇe metody, kter´e se skuteˇcnˇe pouˇz´ıvaj´ı, jsou pomˇernˇe sloˇzit´e k pochopen´ı, takˇze je v tomto z´akladn´ım kurzu nelze dost dobˇre vysvˇetlit. Spokoj´ıme se tedy s konstatov´an´ım, ˇze existuje znaˇcn´e mnoˇzstv´ı v´ykonn´ych iteraˇcn´ıch metod, viz napˇr. [13]. Tak tˇreba pro soustavy s pozitivnˇe definitn´ı matic´ı patˇr´ı mezi nejpopul´arnˇejˇs´ı metoda sdruˇzen´ych gradient˚u (struˇcnˇe CG podle anglick´eho conjugate gradient), v MATLABu 37 viz funkce pcg. Z´akladn´ı verze t´eto metody je struˇcnˇe pops´ana v kapitole 6.2, viz tak´e cviˇcen´ı 6.7. Pro soustavy s nesymetrickou matic´ı je situace komplikovanˇejˇs´ı, jednoznaˇcn´y favorit mezi metodami pro jejich ˇreˇsen´ı neexistuje. Jednou z mnoha pouˇz´ıvan´ych metod je zobecnˇen´a metoda minim´aln´ıch rezidu´ı (struˇcnˇe GMRES podle anglick´eho generalized minimal residual), v MATLABu viz funkce gmres. 0 100 200 300 400 500 600 700 0 200 400 600 800 1000 1200 1400 1600 Rovnic Iterací Jacobi Gauss−Seidel Optimální SOR Obr. 2.3: Srovn´an´ı klasick´ych iteraˇcn´ıch metod 2.3. Cviˇcen´ı 2.1. Metodou GEMz a GEM s ˇc´asteˇcn´ym v´ybˇerem hlavn´ıho prvku (GEMpp) ˇreˇste rovnice a) x+ 2y+ 3z= 6 2x+ 4y+ 5z= 11 7x+ 8y+ 9y= 24 b)10−15x+y= 1 + 10−15 x+ 1011y= 1 + 1011 [ a) GEMz selˇze, GEMpp najde ˇreˇsen´ı (1 1 1)T; b) GEMz spoˇcte (0,888 1)T, GEMpp d´av´a (1 1)T. ] 2.2. Odvod’te (pro n= 3) vzorce (2.13). [LLT=  l11 0 0 l21 l22 0 l31 l32 l33    l11 l21 l31 0l22 l32 0 0 l33  =  l2 11 l11l21 l11l31 l21l11 l2 21 +l2 22 l21l31 +l22l32 l31l11 l31l21 +l32l22 l2 31 +l2 32 +l2 33  , LLT=A⇒l11 =√a11, l21 =a21/l11, l31 =a31/l11, l22 =pa22 −l2 21, . . . ] 2.3. Pomoc´ı Cholesk´eho rozkladu ˇreˇste soustavy rovnic Axi=bi,i= 1,2,3, kde A=  2 1 2 1 2 2 2 2 3 ,b1=  1 −1 0 ,b2=  0 −1 −1 . 38 [L. =  1,414214 0 0 0,707107 1,224745 0 1,414214 0,816497 0,577350 ,x1=  1 −1 0 ,x2=  1 0 −1 .] 2.4. ˇ Reˇste soustavy rovnic Axi=bi,i= 1,2,3, kde A=  1−1 1 2−1 1 1 1 2 ,b1= −1 0 0 ,b2=  1 2 7 ,b3= 2x2−x1. [ Proved’te LU rozklad PA =LU, kde L=  1 0 0 0,5 1 0 0,5−0,33 1 ,U=  2−1 1 0 1,5 1,5 0 0 1 , P=  0 1 0 0 0 1 1 0 0 .Nejdˇr´ıve z´ısk´ame x1,x2, sestav´ıme b3a opˇet LU rozkladem spoˇcteme x3: x1=  1 1 −1 ,x2=  1 2 2 ,b3=  1 3 5 ,x3=  1 1,66 0,66 .] 2.5. Spoˇctˇete ˇc´ıslo podm´ınˇenosti matice soustavy Aze cviˇcen´ı 2.1 b). [κ(A). = 1022 v normˇe k·k∞. ] 2.6. Spoˇctˇete determinant matice Aze cviˇcen´ı 2.4 pouˇzit´ım LU rozkladu. [|A|= (−1)2|U|= 2 ·1,5·1 = 3. ] 2.7. Spoˇctˇete inverzn´ı matici k matici Aze cviˇcen´ı 2.4. [A−1= −1 1 0 −1 0,33 0,33 1−0,66 0,33 . ] 2.8. Naprogramujte metodu GEMz a upravte ji tak, aby bylo moˇzn´e volitelnˇe zapnout/vypnout ˇc´asteˇcn´y v´ybˇer hlavn´ıho prvku (GEMpp). 2.9. Naprogramujte v´ypoˇcet vektorov´ych a maticov´ych norem. 2.10. Naprogramujte ˇreˇsen´ı soustav s tˇr´ıdiagon´aln´ı matic´ı. Prvky matice drˇzte v pamˇeti poˇc´ıtaˇce v ´usporn´em form´atu napˇr. jako tˇri vektory (jednorozmˇern´a pole). Jak vypad´a transformovan´a matice e A:= LU a matice LaUz´ıskan´e LU rozkladem matice A=    100 12 0 0 21 100 23 0 0 32 100 34 0 0 43 100    ? [e A. =    100 12 0 0 0,21 97,48 23 0 0 0,3283 92,4497 34 0 0 0,4651 84,186    , L. =    1 0 0 0 0,21 1 0 0 0 0,3283 1 0 0 0 0,4651 1    ,U. =    100 12 0 0 0 97,48 23 0 0 0 92,4497 34 0 0 0 84,186    . ] 39 2.11. Naprogramujte (a) Jacobiovu a (b) Gauss-Seidelovu iteraˇcn´ı metodu. Jak vypad´a vektor x5pro soustavu rovnic  −194 −4−4 15 33 −3−1   x1 x2 x3 =  29 33 24 , kdyˇz nejdˇr´ıve rovnice soustavy vhodnˇe pˇrerovn´ate (tak, aby matice soustavy byla ryze diagon´alnˇe dominantn´ı) a kdyˇz jako poˇc´ateˇcn´ı aproximaci zvol´ıte vektor x0= (0,0,0)T? [ (a) x5. = (0,995405; 2,007110; 2,986472)T, (b) x5. = (1,000144; 1,999887; 3,000008)T. ] 2.12. Ovˇeˇrte, ˇze pro soustavu rovnic   1 1 1 1 2 2 1 2 3   x1 x2 x3 =  3 5 6  plat´ı: (a) matici soustavy nelze ˇz´adn´ym pˇrerovn´an´ım ˇr´adk˚u upravit tak, aby byla ryze diagon´alnˇe dominantn´ı; (b) Jacobiova metoda nekonverguje (pouˇzijte poˇc´ıtaˇc); (c) matice soustavy je pozitivnˇe definitn´ı (pouˇzijte Sylvesterovo kriterium); (d) Gaussova-Seidelova metoda konverguje (pouˇzijte poˇc´ıtaˇc). 2.13. Ovˇeˇrte, ˇze pro soustavu rovnic   1 1 2 1 2 2 2 2 3   x1 x2 x3 =  4 5 7  plat´ı: (a) matici soustavy nelze ˇz´adn´ym pˇrerovn´an´ım ˇr´adk˚u upravit tak, aby byla ryze diagon´alnˇe dominantn´ı; (b) Jacobiova metoda nekonverguje (pouˇzijte poˇc´ıtaˇc); (c) matice soustavy nen´ı pozitivnˇe definitn´ı (pouˇzijte Sylvesterovo kriterium); (d) Gaussova-Seidelova metoda nekonverguje (pouˇzijte poˇc´ıtaˇc). 2.14. Upravte soustavu ze cviˇcen´ı 2.13 na tvar ATAx =ATbs pozitivnˇe definitn´ı matic´ı soustavy. Jako poˇc´ateˇcn´ı aproximaci zvolte vektor x0= (0,0,0)Ta ˇreˇste Gaussovou-Seidelovou metodou. Pouˇzijte poˇc´ıtaˇc a pozorujte pomalou konvergenci k ˇreˇsen´ı. Kolik iterac´ı je tˇreba k dosaˇzen´ı pˇresnosti ε= 10−6, pouˇzijeme-li k ukonˇcen´ı v´ypoˇctu podm´ınku kxk+1 −xkk∞≤εkxkk∞? Konverguje Jacobiova metoda? [ATA=  6 7 10 7 9 12 10 12 17 ,ATb=  23 28 39 ,k+ 1 = 681. Jacobiova metoda nekonverguje.] 2.15. Naprogramujte Jacobiovu relaxaˇcn´ı metodu (struˇcnˇe JOR) a Gaussovu-Seidelovu relaxaˇcn´ı metodu SOR. Pracujte se soustavami (i)10 3 2 20x1 x2=10 2,(ii)3 2 2 3x1 x2=5 5, poˇc´ateˇcn´ı iteraci zvolte x0= (0,0)T. (a) Jak vypadaj´ı vektory x5vznikaj´ıc´ı pˇri ˇreˇsen´ı soustavy (i) metodou JOR a SOR, kdyˇz pouˇzijeme relaxaˇcn´ı parametr ω= 0,5? (b) Experiment´alnˇe najdˇete pˇribliˇznou optim´aln´ı hodnotu ωopt pro metodu SOR a soustavu (ii). Pro nalezen´e ωopt, a tak´e pro ω= 0,5, urˇcete nejmenˇs´ı poˇcet krok˚u, pˇri kter´em je splnˇena podm´ınka kxk+1 −xkk∞≤10−6kxkk∞. N´avod: pomoc´ı programu mˇeˇnte ωod 1 s krokem 0,01 do 1,99 a urˇcete pˇr´ıpad, kdy vyjde knejmenˇs´ı. [ (a) JOR: x5. = (0,9592; 0,0166)T, SOR: x5. = (0,9640; 0,0083)T; (b) pro ωopt . = 1,15 je k+ 1 = 10, pro ω= 0,5 je k+ 1 = 50. ] 40 3. Aproximace funkc´ı Aproximovat funkci f(x) znamen´a nahradit ji funkc´ı ϕ(x), kter´a je k f(x) v jist´em smyslu bl´ızk´a. P´ıˇseme ϕ(x)≈f(x). Budeme se zab´yvat dvˇema z´akladn´ımi typy aproximace, a to interpolac´ı a metodou nejmenˇs´ıch ˇctverc˚u. Interpolace je takov´a aproximace, pˇri n´ıˇz ϕ(x) nab´yv´a v zadan´ych bodech xipˇredepsan´ych hodnot yi=f(xi). Nˇekdy nav´ıc ˇz´ad´ame, aby funkce ϕafmˇely v bodech xitak´e stejn´e derivace. Interpolaci je vˇenov´an odstavec 3.1. Metoda nejmenˇs´ıch ˇctverc˚u je takov´a aproximace, pˇri n´ıˇz ϕ(x) prokl´ad´ame mezi zadan´ymi body [xi, yi] tak, aby vzd´alenost funkc´ı faϕbyla v jist´em smyslu minim´aln´ı. Je pˇritom charakteristick´e, ˇze funkce ϕbody [xi, yi] neproch´az´ı. Metoda nejmenˇs´ıch ˇctverc˚u je vyloˇzena v odstavci 3.2. Aproximaci ϕ(x) pouˇzijeme k pˇribliˇzn´emu v´ypoˇctu hodnot funkce f(x), tˇreba pˇri vykreslov´an´ı ϕ≈f. Je ˇz´adouc´ı, aby v´ypoˇcet ϕ(x) byl jednoduch´y. Proto se ϕˇcasto hled´a ve tvaru polynomu. Obecnˇe, ϕ(x) se pouˇz´ıv´a k ˇreˇsen´ı ´uloh, v nichˇz vystupuje funkce f, kterou je ´uˇceln´e nebo dokonce nezbytn´e nahradit jej´ı vhodnou aproximac´ı ϕ. Jako pˇr´ıklad uved’me v´ypoˇcet derivace nebo urˇcit´eho integr´alu: f′(x) nahrad´ıme pomoc´ı ϕ′(x) a Rb af(x) dxnahrad´ıme pomoc´ı Rb aϕ(x) dx. 3.1. Interpolace Interpolaˇcn´ı funkci ϕ(x) vyb´ır´ame z vhodn´e tˇr´ıdy funkc´ı. Omez´ıme se na dva nejbˇeˇznˇejˇs´ı pˇr´ıpady: a) ϕ(x) je polynom; b) ϕ(x) je po ˇc´astech polynom, na kaˇzd´em subintervalu obecnˇe jin´y. 3.1.1. Interpolace polynomem Pˇredpokl´adejme, ˇze jsou d´any navz´ajem r˚uzn´e body x0, x1,...,xn, xi6=xjpro i6=j , ˇr´ık´ame jim tak´e uzly interpolace, a v kaˇzd´em z nich je pˇredeps´ana hodnota yi. Hled´ame interpolaˇcn´ı polynom Pn(x) stupnˇe nejv´yˇse n, kter´y splˇnuje interpolaˇcn´ı podm´ınky Pn(xi) = yi, i = 0,1,...,n. (3.1) Existenci interpolaˇcn´ıho polynomu dok´aˇzeme tak, ˇze ho zkonstruujeme. Lagrange˚uv tvar interpolaˇcn´ıho polynomu m´a vyj´adˇren´ı Pn(x) = y0ℓ0(x) + y1ℓ1(x) + ···+ynℓn(x) = n X i=0 yiℓi(x),(3.2) 41 3.1.2. Interpolaˇcn´ı splajny Jestliˇze chceme interpolovat funkci f(x) na pomˇernˇe dlouh´em intervalu ha, bi, mus´ıme ˇz´adat splnˇen´ı interpolaˇcn´ıch podm´ınek v dostateˇcnˇe velk´em poˇctu bod˚u. Pokud bude interpolantem polynom, mus´ı b´yt vysok´eho stupnˇe a to, jak v´ıme, obvykle vede k velk´ym chyb´am mezi uzly. Tudy proto cesta nevede. Lepˇs´ı je rozdˇelit interval ha, bina ˇradu menˇs´ıch subinterval˚u a na kaˇzd´em z nich sestrojit interpolaˇcn´ı polynom niˇzˇs´ıho stupnˇe. Pˇredpokl´adejme, ˇze a=x0< x1<···< xi−1< xi< xi+1 <···< xn−1< xn=b(3.14) je dˇelen´ı intervalu ha, bi. V kaˇzd´em uzlu xije pˇredeps´ana hodnota yiinterpolantu. D´elku i-t´eho intervalu hxi−1, xiioznaˇc´ıme hia d´elku nejdelˇs´ıho intervalu h, tj. hi=xi−xi−1, i = 1,2, . . . , n , h = max 1≤i≤nhi.(3.15) Hledan´y po ˇc´astech polynomick´y interpolant budeme znaˇcit S(x) a nazveme ho interpolaˇcn´ım splajnem. Na kaˇzd´em intervalu hxi−1, xiije S(x) polynom, jehoˇz pˇr´ısluˇsnost ki-t´emu intervalu vyznaˇc´ıme indexem i, tj. S(x) je na intervalu hxi−1, xiipolynom Si(x). K vyj´adˇren´ı polynomu Si(x) s v´yhodou pouˇzijeme lok´aln´ı promˇennou s=x−xi−1. Budeme tak´e pouˇz´ıvat prvn´ı pomˇernou diferenci δi=yi−yi−1 xi−xi−1 =yi−yi−1 hi . Line´arn´ı interpolaˇcn´ı splajn (d´ale jen line´arn´ı splajn) je to nejjednoduˇsˇs´ı, co n´as napadne: kaˇzd´e dva sousedn´ı body [xi−1, yi−1] a [xi, yi] spoj´ıme ´useˇckou. Zˇrejmˇe Si(x) = yi−1+yi−yi−1 xi−xi−1 (x−xi−1) = yi−1+sδi(3.16) je line´arn´ı interpolaˇcn´ı polynom proch´azej´ıc´ı body [xi−1, yi−1] a [xi, yi]. Line´arn´ı splajn S(x) je spojit´a funkce, derivace S′(x) je vˇsak ve vnitˇrn´ıch uzlech obecnˇe nespojit´a. Jestliˇze yi=f(xi), i= 0,1,...,n, a f∈C2ha, bi, pak pro chybu interpolace plat´ı |f(x)−S(x)| ≤ Ch2,(3.17) kde x∈ ha, bije libovoln´e a Cje konstanta nez´avisl´a na h. Pro dostateˇcnˇe mnoho uzl˚u lze uˇcinit chybu libovolnˇe malou. Napˇr´ıklad pˇri vykreslov´an´ı na obrazovku monitoru s rozliˇsen´ım 1920 x 1080 bod˚u jistˇe staˇc´ı pouˇz´ıt 1920 interpolaˇcn´ıch uzl˚u k z´ısk´an´ı kvalitn´ıho grafu interpolovan´e funkce. Moˇzn´a bychom byli spokojeni, i kdybychom zvolili m´enˇe uzl˚u, avˇsak pˇri postupn´em sniˇzov´an´ı poˇctu uzl˚u by nutnˇe nastal okamˇzik, kdy by n´as jiˇz zaˇcaly ruˇsit ostr´e hrany grafu S(x) v interpolaˇcn´ıch 48 uzlech. Pokud bychom souˇcasnˇe vykreslovali tak´e funkci f(x), pak by n´am zaˇcaly vadit tak´e viditeln´e odchylky interpolantu S(x) od interpolovan´e funkce f(x) mezi uzly interpolace. Pˇresnˇejˇs´ı interpolant bychom mohli sestrojit tak, ˇze bychom na intervalech hx0, xki, hxk, x2ki,... aproximovali f(x) pomoc´ı interpolaˇcn´ıch polynom˚u stupnˇe (nejv´yˇse) k > 1. Chyba interpolace by v tom pˇr´ıpadˇe byla ´umˇern´a hk+1, derivace v uzlech xk, x2k,... by vˇsak z˚ustaly nespojit´e. Velk´e kale nem´a smysl pouˇz´ıvat, jinak bychom zase mohli dostat velk´e chyby mezi uzly interpolace a byli bychom zpˇet v situaci, kter´e jsme se pr´avˇe interpolac´ı po ˇc´astech chtˇeli vyhnout. Velmi popul´arn´ı je aproximace po ˇc´astech kubick´ym polynomem, kter´a je nejen spojit´a, ale m´a tak´e spojit´e prvn´ı nebo dokonce i druh´e derivace. Popisu takov´ych aproximac´ı se budeme vˇenovat v n´asleduj´ıc´ıch odstavc´ıch. Hermit˚uv kubick´y interpolaˇcn´ı splajn (d´ale jen Hermit˚uv kubick´y splajn) hled´ame jako funkci S(x), kter´a a) je v intervalu ha, bispojit´a spolu se svou prvn´ı derivac´ı, tj. S∈C1ha, bi, b) splˇnuje interpolaˇcn´ı podm´ınky S(xi) = yi, S′(xi) = di, i = 0,1,...,n, (3.18) kde yi,dijsou pˇredepsan´e funkˇcn´ı hodnoty a derivace, c) je na kaˇzd´em intervalu hxi−1, xii,i= 1,2,...,n, polynom tˇret´ıho stupnˇe. Si(x) je tedy kubick´y Hermit˚uv polynom jednoznaˇcnˇe urˇcen´y podm´ınkami Si(xi−1) = yi−1, S′ i(xi−1) = di−1, Si(xi) = yi, S′ i(xi) = di. Snadno ovˇeˇr´ıme, ˇze tyto podm´ınky jsou splnˇeny pro Si(x) = yi−1+sdi−1+s23δi−2di−1−di hi +s3di−1−2δi+di h2 i .(3.19) Funkce S(x) je spojit´a spolu se svou prvn´ı derivac´ı, druh´a derivace uˇz obecnˇe spojit´a nen´ı. Jestliˇze yi=f(xi), di=f′(xi), i= 0,1,...,n, a f∈C4ha, bi, pak pro chybu interpolace plat´ı |f(x)−S(x)| ≤ Ch4,(3.20) kde x∈ ha, bije libovoln´e a Cje konstanta nez´avisl´a na h. Pokud derivace dinejsou k dispozici, mus´ıme je vypoˇc´ıtat pomoc´ı vhodnˇe zvolen´ych dodateˇcn´ych podm´ınek. Kubick´y interpolaˇcn´ı splajn (d´ale jen kubick´y splajn). Smˇernice dive vnitˇrn´ıch uzlech m˚uˇzeme urˇcit tak, ˇze poˇzadujeme, aby splajn S∈C2ha, bi, tj. aby platilo S′′ i(xi) = S′′ i+1(xi), i = 1,2,...,n−1.(3.21) 49 Derivov´an´ım (3.19) dostaneme S′′ i(x) = (6hi−12s)δi+ (6s−4hi)di−1+ (6s−2hi)di h2 i . Pro x=xije s=hi, takˇze S′′ i(xi) = −6δi+ 2di−1+ 4di hi . Pro x=xi−1je s= 0 a S′′ i(xi−1) = 6δi−4di−1−2di hi . Kdyˇz v posledn´ım vzorci zvˇetˇs´ıme index io jedniˇcku, dostaneme S′′ i+1(xi) = 6δi+1 −4di−2di+1 hi+1 . Dosazen´ım do (3.21) tak dostaneme rovnice hi+1di−1+ 2(hi+1 +hi)di+hidi+1 = 3(hi+1δi+hiδi+1), i = 1,2,...,n−1.(3.22) Jestliˇze pˇredep´ıˇseme okrajov´e podm´ınky S′(a) = da, S′(b) = db,(3.23) pak v soustavˇe (3.22) dosad´ıme v prvn´ı rovnici d0:= daa ˇclen h2dapˇrevedeme na pravou stranu, a v posledn´ı rovnici dosad´ıme dn:= dba ˇclen hn−1dbpˇrevedeme na pravou stranu. Soustavu pak ˇreˇs´ıme a z´ısk´ame zb´yvaj´ıc´ı smˇernice di, i = 1,2,...,n−1. Matice soustavy je tˇr´ıdiagon´aln´ı, diagon´alnˇe dominantn´ı, takˇze soustavu lze snadno vyˇreˇsit GEM upravenou pro soustavy s tˇr´ıdiagon´aln´ı matic´ı. V MATLABu lze pro v´ypoˇcet kubick´eho splajnu pouˇz´ıt funkci spline. Jestliˇze yi=f(xi), i= 0,1,...,n,d0=f′(x0), dn=f′(xn), a kdyˇz f∈C4ha, bi, pak pro chybu interpolace opˇet plat´ı (3.20). Obr´azek 3.3 potvrzuje, ˇze pomoc´ı kubick´eho splajnu lze pro data stejn´a jako v pˇr´ıkladu 3.4 dostat zcela vyhovuj´ıc´ı aproximaci Rungeovy funkce. Kubick´y splajn m´a pozoruhodnou extrem´aln´ı vlastnost, kterou si ted’pop´ıˇseme. Oznaˇc´ıme V={v∈C2ha, bi|v(xi) = yi, i = 0,1,...,n,v′(x0) = d0, v′(xn) = dn} mnoˇzinu vˇsech funkc´ı, kter´e maj´ı v intervalu ha, bispojitou druhou derivaci, proch´azej´ı zadan´ymi body [xi, yi], i= 0,1,...,n, a v krajn´ıch bodech a=x0axn=bjejich derivace nab´yvaj´ı pˇredepsan´ych hodnot d0adn. Pak Rb a[v′′(x)]2dxnab´yv´a na mnoˇzinˇe funkc´ı V sv´e nejmenˇs´ı hodnoty pro kubick´y splajn S(x), tj. plat´ı Zb a [S′′(x)]2dx= min v∈VZb a [v′′(x)]2dx . 50 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1Kubický splajn 1/(1+25x2) Obr. 3.3: Aproximace Rungeovy funkce kubick´ym splajnem Tato vlastnost m´a zaj´ımavou interpretaci v mechanice. Je totiˇz zn´amo, ˇze ohybov´a energie homogenn´ıho izotropn´ıho prutu, jehoˇz stˇrednicov´a ˇc´ara m´a rovnici y=v(x), x∈ ha, bi, m´a pˇribliˇznˇe hodnotu E(v) = cRb a[v′′(x)]2dx, kde cje vhodn´a konstanta. A d´ale je tak´e zn´amo, ˇze prut, kter´y je donucen proch´azet pevn´ymi interpolaˇcn´ımi body [xi, yi], zaujme pozici s minim´aln´ı energi´ı. Extrem´aln´ı vlastnost tedy tvrd´ı, ˇze kubick´y splajn aproximuje stˇrednicovou ˇc´aru takov´eho prutu. Jestliˇze smˇernice daadbv krajn´ıch bodech intervalu ha, binezn´ame, osvˇedˇcil se postup, oznaˇcovan´y v anglicky psan´e literatuˇre jako not a knot. Myˇslenka je jednoduch´a: poˇzadujeme, aby splajn byl jednoduch´ym polynomem tˇret´ıho stupnˇe na prvn´ıch dvou intervalech, tj. pro x0≤x≤x2, a na posledn´ıch dvou intervalech, tj. pro xn−2≤x≤xn. V uzlech x1axn−1tedy uˇz nedoch´az´ı k napojov´an´ı dvou r˚uzn´ych polynom˚u, tj. uzel x1 axn−1uˇz nen´ı knot, ˇcesky uzel splajnu, odtud n´azev postupu not a knot. Polynomy S1(x) a S2(x) maj´ı v bodˇe x1spoleˇcnou funkˇcn´ı hodnotu y1, stejnou prvn´ı derivaci d1a podle (3.21) tak´e stejnou druhou derivaci. Aby oba polynomy byly totoˇzn´e staˇc´ı, kdyˇz budou m´ıt v bodˇe x1tak´e stejnou tˇret´ı derivaci. Stejnou ´uvahu lze prov´est v bodˇe xn−1. Dost´av´ame tak okrajov´e podm´ınky S′′′ 1(x1) = S′′′ 2(x1), S′′′ n−1(xn−1) = S′′′ n(xn−1).(3.24) Kdyˇz pomoc´ı (3.19) vyj´adˇr´ıme podm´ınku S′′′ 1(x1) = S′′′ 2(x1) a uprav´ıme ji pomoc´ı prvn´ı rovnice soustavy (3.22), dostaneme rovnici h2d0+ (h2+h1)d1= [(3h1+ 2h2)h2δ1+h2 1δ2]/(h1+h2).(3.25) Podobnˇe zpracujeme tak´e podm´ınku S′′′ n−1(xn−1) = S′′′ n(xn−1) a dostaneme rovnici (hn+hn−1)dn−1+hn−1dn= [h2 nδn−1+ (2hn−1+ 3hn)hn−1δn]/(hn−1+hn).(3.26) 51 Nezn´am´e d0, d1,...,dnpak dostaneme jako ˇreˇsen´ı soustavy rovnic, z nichˇz prvn´ı je rovnice (3.25), pak n´asleduje n−1 rovnic (3.22) a nakonec pˇrijde jeˇstˇe rovnice (3.26). Pouˇzijeme opˇet GEM upravenou pro soustavy s tˇr´ıdiagon´aln´ı matic´ı. Je-li aproximovan´a funkce periodick´a s periodou b−a, je pˇrirozen´e poˇzadovat, aby aproximuj´ıc´ı splajn byl rovnˇeˇz periodick´y s tout´eˇz periodou. Je-li tedy Speriodick´y kubick´y splajn s periodou b−a, pak na intervalu hb, b +h1inab´yv´a stejn´ych hodnot jako na intervalu ha, a +h1i. Jestliˇze oznaˇc´ıme hn+1 =h1,xn+1 =xn+hn+1,S(xn+1) = yn+1 aS′(xn+1) = dn+1, pak mus´ı platit yn=y0,dn=d0,yn+1 =y1adn+1 =d1. Spojitost S′′ v bodˇe xnznamen´a, ˇze rovnice (3.22) plat´ı tak´e pro i=n. Periodick´y kubick´y splajn je tedy urˇcen, pokud vypoˇcteme d1, d2,...,dnz rovnic (3.22) pro i= 1,2,...,ns t´ım, ˇze v prvn´ı z nich, tj. pro i= 1, m´ısto d0p´ıˇseme dn, v posledn´ı z nich , tj. pro i=n, m´ısto dn+1 p´ıˇseme d1a poloˇz´ıme hn+1 =h1,δn+1 =δ1. Shrnut´ı. Kubick´y interpolaˇcn´ı splajn S(x) je funkce, kter´a a) je v intervalu ha, bispojit´a spolu se svou prvn´ı a druhou derivac´ı, tj. S∈C2ha, bi, b) splˇnuje interpolaˇcn´ı podm´ınky S(xi) = yi,i= 0,1,...,n, kde yijsou pˇredepsan´e funkˇcn´ı hodnoty, c) je na kaˇzd´em intervalu hxi−1, xiipolynom tˇret´ıho stupnˇe, d) splˇnuje okrajov´e podm´ınky (3.23) nebo (3.24) nebo je periodick´a s periodou b−a. Pˇr´ıklad 3.6. Oblouk [x(t), y(t)] ≡[cos t, sin t], t∈ h0, π/2i, budeme aproximovat kˇrivkou [Sx(t), Sy(t)], kde Sx(t) resp. Sy(t) jsou kubick´e splajny funkc´ı cos tresp. sin tna intervalu h0, π/2i. Interval h0, π/2irozdˇel´ıme na nstejn´ych d´ılk˚u, takˇze ti=πi/(2n). Oznaˇc´ıme dx i= [Sx]′(ti), dy i= [Sy]′(ti). Protoˇze x′(t) = −sin t,y′(t) = cos t, poloˇz´ıme dx 0=−sin 0 = 0 , dx n=−sin(π/2) = −1, dy 0= cos 0 = 1 , dy n= cos(π/2) = 0 . V dalˇs´ım budeme uvaˇzovat n= 3. Soustavu rovnic (3.22) sestav´ıme zvl´aˇst ’pro x-ovou a zvl´aˇst ’pro y-ovou sloˇzku. Matice obou soustav je stejn´a, prav´e strany jsou vˇsak r˚uzn´e. Na dalˇs´ım ˇr´adku je uvedena soustava rovnic se dvˇema prav´ymi stranami a jej´ı ˇreˇsen´ı: 1 6π1 4 1 0 0 1 4 1    0 1 dx 1dy 1 dx 2dy 2 −1 0    = 3 −1 2 1 2√3 −1 2√31 2=⇒dy 1=−dx 2 . = 0,865537, dy 2=−dx 1 . = 0,499813. Sx i(t) a Sy i(t) urˇc´ıme podle (3.19). Napˇr´ıklad pro t∈ hπ/6, π/3idostaneme Sx 2(t). = 0,8660 −0,4998(t−1 6π)−0,4431(t−1 6π)2+ 0,1195(t−1 6π)3, Sy 2(t). = 0,5000 + 0,8655(t−1 6π)−0,2554(t−1 6π)2−0,1195(t−1 6π)3. Konstrukce kubick´eho interpolaˇcn´ıho splajnu uˇzit´ım druh´ych derivac´ı. Snadno ovˇeˇr´ıme, ˇze kubick´y polynom Si(x) = yi−1+s6δi−2hiMi−1−hiMi 6+s2Mi−1 2+s3Mi−Mi−1 6hi .(3.27) 52 splˇnuje podm´ınky Si(xi−1) = yi−1, S′′ i(xi−1) = Mi−1, Si(xi) = yi, S′′ i(xi) = Mi. Funkce S(x), kter´a je na kaˇzd´em intervalu hxi−1, xiidefinov´ana pˇredpisem (3.27), proto zˇrejmˇe splˇnuje podm´ınky S(xi) = yi,S′′(xi) = Mi,i= 0,1,...,n.S(x) je tedy na intervalu ha, bispojit´a a m´a v nˇem spojitou druhou derivaci. Abychom dostali kubick´y splajn, mus´ı m´ıt S(x) v ha, bispojitou tak´e prvn´ı derivaci. Protoˇze nespojitost S′(x) m˚uˇze nastat jedinˇe ve vnitˇrn´ıch uzlech, staˇc´ı poˇzadovat S′ i(xi) = S′ i+1(xi), i = 1,2,...,n−1.(3.28) Vyj´adˇr´ıme-li (3.28) pomoc´ı (3.27), dostaneme rovnice hiMi−1+ 2(hi+hi+1)Mi+hi+1Mi+1 = 6(δi+1 −δi), i = 1,2,...,n−1.(3.29) Kdyˇz zvol´ıme okrajov´e podm´ınky S′′(a) = Ma, S′′(b) = Mb,(3.30) dosad´ıme je do (3.29), soustavu rovnic vyˇreˇs´ıme a z´ısk´ame Mi,i= 1,2,...,n−1. Vˇsimnˇete si, ˇze matice soustavy (3.29) je symetrick´a. Je samozˇrejmˇe moˇzn´e uvaˇzovat tak´e jin´e typy okrajov´ych podm´ınek, napˇr. (3.23) nebo (3.24). Kubick´y splajn s vlastnost´ı S′′(a) = S′′(b) = 0 se naz´yv´a pˇrirozen´y kubick´y splajn. Je zn´amo, ˇze pˇrirozen´y kubick´y splajn aproximuje pr˚uhyb prostˇe podepˇren´eho nosn´ıku proch´azej´ıc´ıho body [xi, yi]. Mimaj´ı v´yznam ohybov´ych moment˚u v uzlech [xi, yi]. 3.1.3. Interpolace funkc´ı v´ıce promˇenn´ych Omez´ıme se na pˇr´ıpad, kdy fje funkce dvou promˇenn´ych definovan´a v oblasti Ω. Interpolace po ˇc´astech line´arn´ı. Pˇredpokl´adejme, ˇze Ω je mnoho´uheln´ık. Oblast Ω triangulujeme, tj. vyj´adˇr´ıme ji jako sjednocen´ı troj´uheln´ık˚u T1, T2,...,Tm, z nichˇz kaˇzd´e dva r˚uzn´e bud’to nemaj´ı ˇz´adn´y spoleˇcn´y bod nebo maj´ı spoleˇcn´y vrchol popˇr´ıpadˇe maj´ı spoleˇcnou stranu. Mnoˇzinu T={Tk}m k=1 vˇsech takov´ych troj´uheln´ık˚u naz´yv´ame triangulac´ı oblasti Ω. Vrcholy troj´uheln´ık˚u triangulace oznaˇc´ıme P1= [x1, y1], P2= [x2, y2], ...,Pn= [xn, yn] a nazveme je uzly triangulace. Pˇredpokl´adejme, ˇze v kaˇzd´em uzlu Pije pˇredeps´ana hodnota fi=f(xi, yi) interpolovan´e funkce f. Po ˇc´astech line´arn´ım interpolantem funkce fna oblasti Ω rozum´ıme funkci S, kter´a je v Ω spojit´a, splˇnuje interpolaˇcn´ı podm´ınky S(xi, yi) = fi,i= 1,2,...,n, a kter´a je na kaˇzd´em troj´uheln´ıku Tk∈Tline´arn´ı. Na Tkje tedy z=S(x, y)≡Sk(x, y) rovnice roviny urˇcen´e hodnotami funkce fve vrcholech Tk. Pˇripomeˇnme, ˇze rovnice roviny proch´azej´ıc´ı body [xa, ya, za], [xb, yb, zb] a [xc, yc, zc] m˚uˇze b´yt vyj´adˇrena ve tvaru  x−xay−yaz−za xb−xayb−yazb−za xc−xayc−yazc−za= 0 . 53 Vypoˇc´ıtat hodnotu z=S(x, y) pro (x, y)∈Ω je snadn´e: urˇc´ıme troj´uheln´ık Tk, v nˇemˇz bod [x, y] leˇz´ı, a vypoˇcteme z=Sk(x, y). Chyba interpolace je t´ım menˇs´ı, ˇc´ım jemnˇejˇs´ı triangulaci zvol´ıme. Kdyˇz f∈C2(Ω) (tj. kdyˇz fje v Ω spojit´a spolu se sv´ymi prvn´ımi a druh´ymi parci´aln´ımi derivacemi), pak |f(x, y)−S(x, y)| ≤ Ch2, kde hje nejdelˇs´ı strana troj´uheln´ık˚u triangulace a Cje konstanta nez´avisl´a na h. Interpolace po ˇc´astech biline´arn´ı. Pˇredpokl´adejme, ˇze Ω = ha, bi×hc, dije obd´eln´ık. Pomoc´ı dˇelen´ı a=x0< x1<···< xn=b,c=y0< y1<···< ym=drozloˇz´ıme v´ychoz´ı obd´eln´ık Ω na menˇs´ı obd´eln´ıky Rij ={(x, y)|xi−1≤x≤xi, yj−1≤y≤yj}, i= 1,2,...,n,j= 1,2,...,m. Pˇredpokl´adejme, ˇze v uzlech [xi, yj] jsou pˇredeps´any hodnoty fij =f(xi, yj) funkce f. Na obd´eln´ıku Rij definujeme funkci Sij(x, y) = fi−1,j−1 xi−x xi−xi−1 yj−y yj−yj−1 +fi,j−1 x−xi−1 xi−xi−1 yj−y yj−yj−1 + +fi−1,j xi−x xi−xi−1 y−yj−1 yj−yj−1 +fij x−xi−1 xi−xi−1 y−yj−1 yj−yj−1 . Funkce Sij je biline´arn´ı , tj. pro pevn´e x=Cje Sij(C, y) line´arn´ı funkce promˇenn´e ya pro pevn´e y=Dje Sij(x, D) line´arn´ı funkce promˇenn´e x. Funkce Sij je interpolant funkce fna obd´eln´ıku Rij , tj. Sij nab´yv´a ve vrcholech [xi−1, yj−1], [xi, yj−1], [xi, yj] a [xi−1, yj] stejn´ych hodnot jako funkce f, jak se snadno pˇresvˇedˇc´ıme. Na Ω definujeme funkci Spˇredpisem S(x, y) = Sij(x, y) pro (x, y)∈Rij,i= 1,2,...,n, j= 1,2,...,m. Protoˇze S(xi, yj) = fij,i= 0,1,...,n,j= 0,1,...,m, ˇrekneme, ˇze Sje po ˇc´astech biline´arn´ı interpolant funkce fna obd´eln´ıku Ω. Snadno ovˇeˇr´ıme, ˇze Sje v Ω spojit´a (staˇc´ı si uvˇedomit, ˇze Sij, a tedy tak´e S, je na kaˇzd´e stranˇe obd´eln´ıka Rij line´arn´ı funkce jednoznaˇcnˇe urˇcena pomoc´ı hodnot funkce fv koncov´ych bodech t´eto strany). Kdyˇz f∈C2(Ω), pak pro chybu interpolace plat´ı odhad |S(x, y)−f(x, y)| ≤ C(h2+k2),kde h= max 1≤i≤n(xi−xi−1), k = max 1≤j≤m(yj−yj−1) aCje konstanta nez´avisl´a na h,k. 3.2. Metoda nejmenˇs´ıch ˇctverc˚u oznaˇcuje postup pro pˇribliˇzn´e ˇreˇsen´ı pˇreurˇcen´ych nebo nepˇresnˇe zadan´ych soustav rovnic, zaloˇzen´y na minimalizaci kvadr´at˚u jejich rezidu´ı. Prokl´ad´an´ı dat kˇrivkami je v´yznamn´a skupina ´uloh, kter´e lze metodou nejmenˇs´ıch ˇctverc˚u ˇreˇsit (v anglicky psan´e literatuˇre se pro tyto aplikace pouˇz´ıv´a oznaˇcen´ı curve fitting). Popiˇsme si, o co v takov´ych ´uloh´ach jde. 54 Necht ’tje nez´avisle promˇenn´a, napˇr´ıklad ˇcas, a y(t) je nezn´am´a funkce promˇenn´e t, kterou chceme aproximovat. Pˇredpokl´adejme, ˇze jsme provedli mpozorov´an´ı, pˇri nichˇz byly hodnoty ypˇribliˇznˇe zmˇeˇreny pro urˇcit´e (navz´ajem r˚uzn´e) hodnoty t, takˇze yi≈y(ti), i = 1,2,...,m, kde symbol ≈vyjadˇruje pˇribliˇznou rovnost. Naˇs´ım z´amˇerem je modelovat y(t) line´arn´ı kombinac´ı nb´azov´ych funkc´ı pro nˇejak´e n≤m: y(t)≈x1ϕ1(t) + x2ϕ2(t) + ···+xnϕn(t) =: Rn(t). Funkce Rn(t) se ve statistice naz´yv´a line´arn´ı regresn´ı funkce. B´azov´e funkce navrhujeme podle oˇcek´avan´eho pr˚ubˇehu nezn´am´e funkce y(t), urˇcit se maj´ı parametry x1, x2,...,xn, a to tak, aby yi≈Rn(ti), i = 1,2,...,m, maticovˇe y≈Ax, kde y= (y1, y2,...,ym)Tjsou namˇeˇren´a data, x= (x1, x2,...,xn)Tje vektor nezn´am´ych parametr˚u a Aje tak zvan´a n´avrhov´a matice, A=       ϕ1(t1)ϕ2(t1). . . ϕn(t1) ϕ1(t2)ϕ2(t2). . . ϕn(t2) . . .. . .. . . ϕ1(tm)ϕ2(tm). . . ϕn(tm)       ≡(ϕ1,ϕ2,...,ϕn). Vektor ϕi= (ϕi(t1), ϕi(t2),...ϕi(tm))T,i= 1,2,...,n, je i-t´y sloupec matice A. Rezidua jsou rozd´ıly mezi pozorov´an´ımi yia modelovan´ymi hodnotami Rn(ti): ri=yi−Rn(ti) = yi− n X j=1 ϕj(ti)xj≡yi− n X j=1 aijxj, i = 1,2,...,m, kde aij =ϕj(ti). V maticov´em z´apisu r=y−Ax .(3.31) Parametry xichceme urˇcit tak, aby rezidua byla co nejmenˇs´ı. Metodu nejmenˇs´ıch ˇctverc˚u dostaneme, kdyˇz minimalizujeme souˇcet ˇctverc˚u rezidu´ı: krk2 2= m X i=1 r2 i→min .(3.32) Nˇekdy se pouˇz´ıv´a tak´e v´aˇzen´a metoda nejmenˇs´ıch ˇctverc˚u: kdyˇz jsou nˇekter´a pozorov´an´ı v´yznamnˇejˇs´ı neˇz ostatn´ı, m˚uˇzeme jednotliv´ym pozorov´an´ım pˇrisoudit v´ahy wi>0 a minimalizovat souˇcet v´aˇzen´ych ˇctverc˚u rezidu´ı krk2 2,w = m X i=1 wir2 i→min . 55 ri=yi−Rn (ti ) r1 2 r2 2 ri−1 2 ri 2 ri+1 2 rm 2 Σi=1 m ri 2 → min t1 y1 t2 y2 ti−1 yi−1 ti yi ti+1 yi+1 tm ym Obr. 3.4: Princip metody nejmenˇs´ıch ˇctverc˚u Je-li napˇr´ıklad chyba i-t´eho pozorov´an´ı pˇribliˇznˇe rovna ei, zvol´ıme wi= 1/ei. Norm´aln´ı soustava rovnic. Oznaˇcme F(x) = ky−Axk2 2≡ krk2 2.ˇ Reˇsen´ı minimalizaˇcn´ı ´ulohy (3.32) mus´ı splˇnovat nutnou podm´ınku pro extr´em: ∂F(x) ∂xk =∂ ∂xk m X i=1 yi− n X j=1 aijxj!2 = 0 , k = 1,2,...,n. Kdyˇz provedeme naznaˇcen´e derivov´an´ı, dostaneme ∂F(x) ∂xk = 2 m X i=1 yi− n X j=1 aijxj!(−aik) = 0, a odtud n X j=1 m X i=1 aikaij!xj= m X i=1 aikyi, k = 1,2,...,n, coˇz lze zapsat maticovˇe jako ATAx =ATy.(3.33) Soustava line´arn´ıch rovnic (3.33) je zn´ama jako norm´aln´ı soustava rovnic. Kdyˇz jsou sloupce matice Aline´arnˇe nez´avisl´e, je matice G:= ATApozitivnˇe definitn´ı, takˇze ˇreˇsen´ı x∗norm´aln´ı soustavy rovnic minimalizuje F(x)≡ krk2 2a je tedy ˇreˇsen´ım ´ulohy (3.32): ky−Ax∗k2 2= min x∈Rnky−Axk2 2nebo-li x∗= argmin x∈Rnky−Axk2 2. 56 Vyj´adˇr´ıme-li norm´aln´ı soustavu rovnic pomoc´ı vektor˚u ϕi, dostaneme       (ϕ1,ϕ1) (ϕ1,ϕ2)... (ϕ1,ϕn) (ϕ2,ϕ1) (ϕ2,ϕ2)... (ϕ2,ϕn) . . .. . .... . . . (ϕn,ϕ1) (ϕn,ϕ2) (ϕn,ϕn)             x1 x2 . . . xn       =      (ϕ1,y) (ϕ2,y) . . . (ϕn,y)       ,(3.34) kde (ϕk,ϕj) = m X i=1 ϕk(ti)ϕj(ti) a (ϕk,y) = m X i=1 ϕk(ti)yi jsou skal´arn´ı souˇciny vektor˚u ϕk,ϕjaϕk,y. Matice Gsoustavy (3.34) se naz´yv´a Gramova matice soustavy vektor˚u {ϕj}n j=1. Pˇri n´avrhu aproximace Rn(t) bychom mˇeli vyb´ırat funkce ϕi(t) tak, aby sloupce ϕi matice Abyly line´arnˇe nez´avisl´e. V opaˇcn´em pˇr´ıpadˇe, jak se d´a uk´azat, m´a ´uloha (3.32) nekoneˇcnˇe mnoho ˇreˇsen´ı, coˇz je zˇrejmˇe neˇz´adouc´ı. Uved’me si dva v´yznamn´e speci´aln´ı pˇr´ıpady, pro kter´e jsou sloupce matice Aline´arnˇe nez´avisl´e (d˚ukaz lze naj´ıt napˇr. v [1]): a) pro n=N+ 1 vol´ıme ϕj(t) = tj−1,j= 1,2,...,N + 1; b) pro n= 2N+ 1 vol´ıme ϕ1(t) = 1, ϕ2k(t) = cos kπ Lt, ϕ2k+1(t) = sin kπ Lt, k = 1,2,...,N, a ˇcasy pozorov´an´ı tivyb´ır´ame z intervaluhc, c + 2L), kde L > 0, clibovoln´e. Aproximace Rn(t) je v pˇr´ıpadˇe a) algebraick´y polynom stupnˇe Na v pˇr´ıpadˇe b) trigonometrick´y polynom stupnˇe N. Kdyˇz je m=na matice Aje regul´arn´ı, pak x∗=A−1yar=o, tj. Rn(ti) = yi, i= 1,2,...,m. Pokud jsou vˇsak namˇeˇren´a data yizat´ıˇzena chybami, pak nen´ı ´uˇceln´e, aby funkce Rn(t) tyto chyby kop´ırovala. Naopak, aby Rn(t) vˇerohodnˇe vystihovala (rekonstruovala) nezn´amou funkci y(t), je ˇz´adouc´ı, aby Rn(t) namˇeˇren´a data vyrovn´avala (vyhlazovala). To je ale moˇzn´e jen tehdy, kdyˇz poˇcet pozorov´an´ı mje v´yraznˇe vˇetˇs´ı neˇz poˇcet nn´avrhov´ych parametr˚u, tj. pro m≫n. Pˇr´ıklad 3.7. Pro data pˇredepsan´a tabulkou ti0 0,5 1 1,5 2 2,5 3 yi3,57 2,99 2,62 2,33 2,22 2,10 2,05 urˇc´ıme aproximaci R2(t) = x1+x2e−tmetodou nejmenˇs´ıch ˇctverc˚u. Zˇrejmˇe ϕ1(t) = 1 aϕ2(t) = e−t. Norm´aln´ı soustava rovnic je tvaru     7 P i=1 1·1 7 P i=1 1·e−ti 7 P i=1 e−ti·1 7 P i=1 e−ti·e−ti        x1 x2    =    7 P i=1 1·yi 7 P i=1 e−ti·yi    . 57 f′(x) = f(x+h)−f(x−h) 2h−1 6h2f′′′(ξ),(4.3b) f′(x) = 3f(x)−4f(x−h) + f(x−2h) 2h+1 3h2f′′′(ξ).(4.3c) Uved’me jeˇstˇe nejzn´amˇejˇs´ı formuli f′′(x1). =P′′ 2(x1) pro v´ypoˇcet druh´e derivace. Rovnost f′′(x) = f(x+h)−2f(x) + f(x−h) h2−1 12h2f(4)(ξ).(4.4) ovˇeˇr´ıme uˇzit´ım Taylorova rozvoje f(x±h) okolo x. Formule ze vzorce (4.2a) je zn´ama jako prvn´ı diference vpˇred (dopˇredn´a diference) a formule ze vzorce (4.2b) jako prvn´ı diference vzad (zpˇetn´a diference). Formule ze vzorce (4.3b) b´yv´a oznaˇcov´ana jako prvn´ı centr´aln´ı diference a formule ze vzorce (4.4) jako druh´a centr´aln´ı diference. Numerick´y v´ypoˇcet parci´aln´ı derivace nepˇredstavuje ˇz´adn´y nov´y probl´em: derivujeme-li podle promˇenn´e xi, ostatn´ıch promˇenn´ych xj6=xisi nevˇs´ım´ame a nˇekterou z v´yˇse uveden´ych formul´ı aplikujeme jen na xi. Tak tˇreba pomoc´ı dopˇredn´e diference (4.2a) dostaneme ∂f(x1, x2) ∂x2≈f(x1, x2+h)−f(x1, x2) h. Podm´ınˇenost numerick´eho v´ypoˇctu derivace. Ve vzorc´ıch (4.2) – (4.4) jsme uvedli vˇzdy formuli (jako prvn´ı sˇc´ıtanec na prav´e stranˇe) a jej´ı diskretizaˇcn´ı chybu (jako druh´y sˇc´ıtanec). Pˇri numerick´em v´ypoˇctu derivace hraj´ı v´yznamnou roli tak´e zaokrouhlovac´ı chyby, a to jak v hodnot´ach funkce f(tj. ve vstupn´ıch datech), tak pˇri vyhodnocen´ı formule (tj. pˇri v´ypoˇctu). Uk´aˇzeme si to pro formuli ze vzorce (4.2a). Ve skuteˇcnosti za pˇribliˇznou hodnotu derivace f′(x) povaˇzujeme v´yraz ˜ f′(x) := ˜ f(x+h)−˜ f(x) h=[f(x+h) + ε1]−[f(x) + ε0] h=f′(x)+ 1 2hf′′(ξ)+ ε1−ε0 h, kde ε1resp. ε0je zaokrouhlovac´ı chyba, kter´e se dopust´ıme pˇri v´ypoˇctu f(x+h) resp. f(x). Tedy f′(x) = ˜ f(x+h)−˜ f(x) h+Ed+Er, kde Ed:= −1 2hf′′(ξ) je diskretizaˇcn´ı chyba aEr:= −(ε1−ε0)/h je chyba zaokrouhlovac´ı. Chov´an´ı obou chyb je pro h→0 diametr´alnˇe odliˇsn´e: zat´ımco |Ed| → 0, |Er| → ∞. Pro mal´a hse tedy zˇrejmˇe jedn´a o ˇspatnˇe podm´ınˇenou ´ulohu: mal´e zmˇeny ε0, ε1ve vstupn´ıch datech vyvolaj´ı velkou zmˇenu Era n´aslednˇe tak´e v´ysledku ˜ f′(x). Kdyˇz pro jednoduchost zanedb´ame zaokrouhlovac´ı chyby vznikaj´ıc´ı pˇri vyˇc´ıslen´ı formule [ ˜ f(x+h)−˜ f(x)]/h, dost´av´ame pro celkovou chybu E=Ed+Erodhad |E| ≤ |Ed|+|Er| ≤ 1 2hM2+ 2 ε h≡g(h), 64 kde M2≥ |f′′(ξ)|aε≥max(|ε1|,|ε2)|). Minimalizac´ı funkce g(h) obdrˇz´ıme optim´aln´ı d´elku kroku hopt = 2rε M2 ,pro kterou |Eopt|=g(hopt) = 2pεM2. Pˇredpokl´adejme, ˇze hodnoty f(x) i f(x+h) dok´aˇzeme vypoˇc´ıtat s relativn´ı chybou rovnou pˇribliˇznˇe ˇc´ıslu δ, takˇze ε≈M0δ, kde M0≈max(|f(x0)|,|f(x0+h)|). Pro M0≈M2 je hopt ≈2√δa|Eopt| ≈ 2M0√δ. Poˇc´ıt´ame-li tedy napˇr. ve dvojn´asobn´e pˇresnosti a pokud δ≈10−16, pak hopt ≈2·10−8. Jestliˇze nav´ıc |f′(x)| ≈ M0, pak |Eopt| ≈ 2|f′(x)|√δ, a to znamen´a, ˇze velikost relativn´ı chyby derivace ˜ f′(x) je ˇr´adovˇe rovna druh´e odmocninˇe velikosti relativn´ı chyby funkˇcn´ıch hodnot. To n´as opravˇnuje k tvrzen´ı: pˇri pˇribliˇzn´em v´ypoˇctu derivace formul´ı dopˇredn´e (nebo zpˇetn´e) diference doch´az´ı pˇri optim´aln´ı volbˇe kroku ke ztr´atˇe pˇribliˇznˇe poloviny platn´ych cifer. 0 0.2 0.5 1 1.5 2 x 10−7 0 0.5 1 1.5 2x 10−7 h g(h) [hopt,Eopt] Obr. 4.1: Chyba numerick´e derivace: pro g(h) = 1 2h+ 2 ·10−16/h je hopt = 2 ·10−8=Eopt Podobn´e chov´an´ı vykazuj´ı i ostatn´ı formule numerick´eho derivov´an´ı, tj. pro krok h bl´ızk´y hopt je numerick´y v´ypoˇcet derivace ˇspatnˇe podm´ınˇen´a ´uloha: nepatrn´e zmenˇsen´ı kroku vyvol´a znaˇcn´y n´arust chyby, viz obr. 4.1. 4.2. Richardsonova extrapolace Pˇribliˇzn´y v´ypoˇcet derivace lze efektivnˇe zpˇresnit technikou zn´amou jako Richardsonova extrapolace. Je to univerz´aln´ı postup umoˇzˇnuj´ıc´ı pomoc´ı z´akladn´ı metody niˇzˇs´ı pˇresnosti vytv´aˇret metody vyˇsˇs´ı pˇresnosti. Ukaˇzme si, jak se to dˇel´a. Pˇredpokl´adejme, ˇze z´akladn´ı metoda je reprezentov´ana funkc´ı F(h) parametru h. Metodou Fum´ıme vypoˇc´ıtat hodnotu F(h) pro mal´a h > 0. Naˇs´ım c´ılem je co nejpˇresnˇeji aproximovat hodnotu F(0), kterou vˇsak pˇr´ımo z formule Furˇcit neum´ıme. 65 Pˇredpokl´adejme, ˇze funkce F(h) m˚uˇze b´yt zaps´ana ve tvaru mocninn´eho rozvoje F(h) = a0+a1h2+a2h4+a3h6+... (4.5) Konkr´etn´ı formule Fje zkoum´ana v pˇr´ıkladu 4.1, viz (4.14), v pˇr´ıkladu 4.5 a ve cviˇcen´ı 4.2, viz (4.30). Pro mal´e hje F(h) jistˇe dobrou aproximac´ı F(0) = a0. Pokus´ıme se naj´ıt lepˇs´ı aproximaci a0. Zaˇcneme t´ım, ˇze vypoˇcteme F(h 2). Podle (4.5) plat´ı Fh 2=a0+a1h 22 +a2h 24 +a3h 26 +... (4.6) Nejvˇetˇs´ı chybu ve v´yrazu a0−F(h) i a0−F(h 2) pˇredstavuje ˇclen obsahuj´ıc´ı druhou mocninu h. Zbav´ıme se ho tak, ˇze od ˇctyˇrn´asobku rovnice (4.6) odeˇcteme rovnici (4.5) a v´ysledek dˇel´ıme tˇremi. Tak dostaneme F2(h) := 4F(h 2)−F(h) 3=a0+a(2) 2h4+a(2) 3h6+... (4.7) Snadno ovˇeˇr´ıme, ˇze |a(2) j|<|aj|,j= 2,3,....F2(h) je proto lepˇs´ı aproximace a0neˇz F(h) nebot ’a0−F2(h) zaˇc´ın´a aˇz ˇctvrtou mocninou h. Dostali jsme tedy metodu F2, kter´a je (pro dosti mal´a h) lepˇs´ı neˇz p˚uvodn´ı metoda F. Protoˇze F2(h)≈F(0) je spoˇctena pomoc´ı hodnot funkce Fpro hah 2, pˇredstavuje F2(h)extrapolaci funkce Fdo nuly (ovˇeˇrte, ˇze F2(h) = P1(0), kde P1(t) je line´arn´ı interpolaˇcn´ı polynom proch´azej´ıc´ı body [(h 2)2, F(h 2)] a [h2, F(h)]). Podobn´ym postupem odstran´ıme z F2(h) ˇclen obsahuj´ıc´ı ˇctvrtou mocninu ha z´ısk´ame jeˇstˇe lepˇs´ı aproximaci F(0). Nejprve vypoˇcteme F2(h 2). Podle (4.7) plat´ı F2h 2=a0+a(2) 2h 24 +a(2) 3h 26 +... (4.8) Rovnici (4.8) n´asob´ıme 16, odeˇcteme (4.7) a v´ysledek dˇel´ıme 15. Tak dostaneme metodu F3, kter´a je pro zvolen´e hdefinov´ana pˇredpisem F3(h) := 16F2(h 2)−F2(h) 15 =a0+a(3) 3h6+... (4.9) pˇriˇcemˇz |a(3) j|<|a(2) j|<|aj|,j= 3,4,.... Vˇsimnˇete si, abychom mohli vypoˇc´ıtat F2(h 2), mus´ıme nejdˇr´ıve urˇcit F(h 4). Takto m˚uˇzeme pokraˇcovat a z´ısk´avat st´ale lepˇs´ı metody, pro kter´e Fi+1(h) = 4iFi(h 2)−Fi(h) 4i−1=a0+a(i+1) i+1 h2i+2 +. . . , i = 1,2,... (4.10) a kde F1(h) = F(h). Pro koeficienty rozvoje pˇritom plat´ı |a(i+1) j|<|aj|,j=i+1, i+2,.... Protoˇze Fi(h)−F(0) = a(i) ih2i+..., ˇrekneme, ˇze Fi(h) je aproximace F(0) ˇr´adu h2i. 66 V´ypoˇcet lze pˇrehlednˇe uspoˇr´adat do tabulky F1(h) F1(h 2)F2(h) F1(h 4)F2(h 2)F3(h) F1(h 8)F2(h 4)F3(h 2)F4(h) F1(h 16 )F2(h 8)F3(h 4)F4(h 2)F5(h) . . .. . .. . .. . .. . .... ⇐⇒ T00 T10 T11 T20 T21 T22 T30 T31 T32 T33 T40 T41 T42 T43 T44 . . .. . .. . .. . .. . .... Tab. 4.1 Richardsonova extrapolace, Tsi =Fi+1(h/2s−i) Tabulku vyplˇnujeme po ˇr´adc´ıch. Prvky tabulky oznaˇc´ıme jako Tsi, kde s= 0,1,... je ˇr´adkov´y index a i= 0,1,...s je index sloupcov´y. Prvek Ts0v prvn´ım sloupci tabulky vypoˇcteme pomoc´ı z´akladn´ı metody F=F1, Ts0=F(h/2s), s = 0,1,... , (4.11) a dalˇs´ı prvky v tomto ˇr´adku poˇc´ıt´ame ve shodˇe s (4.10) podle pˇredpisu Tsi := 4iTs,i−1−Ts−1,i−1 4i−1=Ts,i−1+Ts,i−1−Ts−1,i−1 4i−1, i = 1,2,...,s. (4.12) V´ypoˇcet ukonˇc´ıme a Tsi povaˇzujeme za dostateˇcnˇe pˇresnou aproximaci F(0), pokud |Tsi −Ts,i−1|<max(εr|Tsi|, εa),(4.13) kde εrje poˇzadovan´a relativn´ı pˇresnost a εapoˇzadovan´a pˇresnost absolutn´ı. Pˇr´ıklad 4.1. Richardsonovu extrapolaci pouˇzijeme pro zpˇresnˇen´ı v´ypoˇctu derivace podle formule (4.3b). Jestliˇze m´a funkce fdostateˇcn´y poˇcet spojit´ych derivac´ı, pak F(h) := f(x+h)−f(x−h) 2h=f′(x) + f(3)(x) 3! h2+f(5)(x) 5! h4+. . . , (4.14) takˇze F(h) je tvaru (4.5). Poˇc´ıtejme derivaci funkce f(x) = cos xpro x= 1. Zvol´ıme napˇr. h= 0,8 a v´ypoˇcet ukonˇc´ıme, kdyˇz bude splnˇena podm´ınka (4.13) pro εr=εa= 10−5. V n´asleduj´ıc´ı tabulce znaˇc´ıme hs=h/2s, prvky Ts0poˇc´ıt´ame ze vztahu Ts0=cos(1 + hs)−cos(1 −hs) 2hs , prvky Ts1aTs2v dalˇs´ıch sloupc´ıch poˇc´ıt´ame podle (4.12). ˇ C´ısla v tabulce jsou zaokrouhlena na 6 cifer. Protoˇze |T32 −T31|<10−5, povaˇzujeme T32 =−0,841471 za pˇribliˇznou hodnotu f′(1). Pˇresn´a hodnota f′(1) = −sin(1) . =−0,84147098, takˇze T32 m´a vˇsechny 67 s hsTs0Ts1Ts2 0 0,8−0,754543 1 0,4−0,819211 −0,840766 2 0,2−0,835872 −0,841426 −0,841470 3 0,1−0,840069 −0,841468 −0,841471 cifry platn´e.  Pozn´amka. Richardsonovu extrapolaci lze aplikovat na z´akladn´ı metodu Ftak´e v pˇr´ıpadˇe, kdyˇz m´a funkce F(h) obecn´y rozvoj F(h) = a0+a1hp1+a2hp2+a3hp3+... (4.5’) kde 1 ≤p1< p2< p3< . . . jsou pˇrirozen´a ˇc´ısla. Pˇresnˇejˇs´ı metodu Fi+1 v tom pˇr´ıpadˇe definujeme pˇredpisem Fi+1(h) = 2piFi(h 2)−Fi(h) 2pi−1=a0+a(i+1) i+1 hpi+1 +. . . , i = 1,2,... , (4.10’) aTsi poˇc´ıt´ame podle Tsi := 2piTs,i−1−Ts−1,i−1 2pi−1=Ts,i−1+Ts,i−1−Ts−1,i−1 2pi−1, i = 1,2,...,s. (4.12’) Protoˇze Fi(h)−F(0) = a(i) ihpi+..., ˇrekneme, ˇzeFi(h) je aproximace F(0) ˇr´adu hpi. Pro pi= 2idostaneme dˇr´ıve uvaˇzovan´y pˇr´ıpad, viz (4.5), (4.10) a (4.12).  Pˇr´ıklad 4.2. Richardsonovou extrapolac´ı zpˇresn´ıme v´ypoˇcet derivace podle formule (4.2a). Z Taylorovy vˇety dostaneme F(h) := f(x+h)−f(x) h=f′(x) + f(2)(x) 2! h+f(3)(x) 3! h2+. . . , (4.15) coˇz odpov´ıd´a (4.5’) pro pi=i. Poˇc´ıtat budeme stejnou ´ulohu jako v pˇr´ıkladu 4.1. Tentokr´at poˇzadovanou pˇresnost dos´ahneme aˇz pro T44. Richardsonova extrapolace je m´enˇe s hsTs0Ts1Ts2Ts3Ts4 0 0,8−0,959381 1 0,4−0,925838 −0,892295 2 0,2−0,889723 −0,853608 −0,840712 3 0,1−0,867062 −0,844401 −0,841332 −0,841421 4 0,05 −0,854625 −0,842188 −0,841451 −0,841468 −0,841471 ´uˇcinn´a: zat´ımco pro formuli (4.3b) je T32 aproximace ˇr´adu h6, pro formuli (4.2a) je T44 aproximace ˇr´adu h5a k dosaˇzen´ı poˇzadovan´e pˇresnosti bylo tˇreba zvolit menˇs´ı hs. 68 4.3. Numerick´e integrov´an´ı C´ılem tohoto odstavce je pˇribliˇzn´y v´ypoˇcet integr´alu I(f) := Rb af(x) dx. Existuje nˇekolik d˚uvod˚u, proˇc tento integr´al nepoˇc´ıt´ame pˇresnˇe, napˇr´ıklad a) integr´al I(f) neum´ıme spoˇc´ıtat analytick´ymi metodami; b) analytick´y v´ypoˇcet je pˇr´ıliˇs pracn´y; c) funkce f(x) je d´ana jen tabulkou. Za pˇribliˇznou hodnotu integr´alu I(f) povaˇzujeme integr´al Q(f) := I(ϕ), kde ϕ(x) je vhodn´a aproximace funkce f(x). Pˇredpis Q(f) pro pˇribliˇzn´y v´ypoˇcet integr´alu se naz´yv´a kvadraturn´ı formule. Rozd´ıl I(f)−Q(f) oznaˇc´ıme R(f) a nazveme (diskretizaˇcn´ı) chybou kvadraturn´ı formule, tedy I(f) = Q(f) + R(f). ˇ Rekneme, ˇze kvadraturn´ı formule je ˇr´adu r, kdyˇz integruje pˇresnˇe polynomy stupnˇe r a polynomy stupnˇe r+ 1 uˇz pˇresnˇe neintegruje, tj. kdyˇz R(xk) = 0 pro k= 0,1,...,r, aR(xr+1)6= 0. ˇ R´ad formule staˇc´ı ovˇeˇrit na intervalu ha, bi=h0,1i. Skuteˇcnˇe, pomoc´ı transformace x=a+t(b−a) dostaneme Rb axkdx= (b−a)R1 0g(t) dt, kde g(t) = (a+t(b−a))kje polynom stupnˇe kv promˇenn´e t. Proto kdyˇz formule integruje pˇresnˇe polynomy stupnˇe k≤rna intervalu h0,1i, integruje pˇresnˇe tak´e polynomy stupnˇe rna intervalu ha, bi. 4.3.1. Z´akladn´ı formule dostaneme integrac´ı interpolaˇcn´ıho polynomu. Obd´eln´ıkov´a formule. Kdyˇz P0(x) = f(1 2(a+b)) je polynom stupnˇe 0, tj. konstanta rovn´a hodnotˇe funkce fve stˇredu 1 2(a+b) intervalu ha, bi, pak odpov´ıdaj´ıc´ı formule QM(f) := Zb a P0(x) dx= (b−a)fa+b 2.(4.16) N´azev formule vyjadˇruje skuteˇcnost, ˇze pro f(1 2(a+b)) >0 je QM(f) obsah obd´eln´ıka o stran´ach d´elky b−aaf(1 2(a+b)). Obd´eln´ıkov´a formule (4.16) se v anglicky psan´e literatuˇre oznaˇcuje jako midpoint rule, odtud index M. Obd´eln´ıkov´a formule je ˇr´adu 1: na intervalu h0,1ije QM(1) = 1 = I(1), QM(x) = 1 2=I(x) a QM(x2) = 1 46=1 3=I(x2). Pro chybu RM(f) obd´eln´ıkov´e formule lze za pˇredpokladu, ˇze f∈C2ha, bi, odvodit RM(f) = 1 24f′′(ξ)(b−a)3,kde ξ∈(a, b) (4.17) je nˇejak´y (bl´ıˇze neurˇcen´y) bod intervalu (a, b). Lichobˇeˇzn´ıkov´a formule. Jako P1(x) oznaˇc´ıme line´arn´ı polynom proch´azej´ıc´ı body [a, f(a)] a [b, f(b)]. Integrac´ı P1(x) na intervalu ha, biobdrˇz´ıme QT(f) := Zb a P1(x) dx=b−a 2[f(a) + f(b)] .(4.18) 69 a (a+b)/2 b a b a (a+b)/2 b Obr. 4.2: Obd´eln´ıkov´a, lichobˇeˇzn´ıkov´a a Simpsonova formule N´azev formule vyjadˇruje skuteˇcnost, ˇze pro f(a)>0, f(b)>0 je QT(f) obsah lichobˇeˇzn´ıka, jehoˇz rovnobˇeˇzn´e strany maj´ı d´elky f(a), f(b) a jehoˇz v´yˇska je rovna b−a. Index Tje prvn´ı p´ısmeno anglick´eho sl˚uvka trapezoid, ˇcesky lichobˇeˇzn´ık. Lichobˇeˇzn´ıkov´a formule je ˇr´adu 1: line´arn´ı polynom integruje pˇresnˇe, kvadratick´y nikoliv. Kdyˇz f∈C2ha, bi, pak pro chybu lichobˇeˇzn´ıkov´e formule plat´ı RT(f) = −1 12f′′(ξ)(b−a)3,kde ξ∈(a, b).(4.19) Vˇsimnˇete si: a) Pokud se druh´a derivace f′′(x) funkce f(x) na intervalu ha, bipˇr´ıliˇs nemˇen´ı, pak je absolutn´ı hodnota |RT(f)|chyby lichobˇeˇzn´ıkov´e fomule pˇribliˇznˇe dvakr´at vˇetˇs´ı neˇz absolutn´ı hodnota |RM(f)|chyby formule obd´eln´ıkov´e. b) Pokud druh´a derivace f′′(x) funkce f(x) nemˇen´ı na intervalu ha, biznam´enko, tj. jeli funkce f(x) poˇr´ad konvexn´ı nebo konk´avn´ı, pak znam´enko chyby lichobˇeˇzn´ıkov´e formule je opaˇcn´e neˇz znam´enko chyby formule obd´eln´ıkov´e. Za tˇechto okolnost´ı pˇresn´a hodnota I(f) integr´alu leˇz´ı v intervalu, jehoˇz krajn´ı body jsou hodnoty QM(f) a QT(f). Simpsonova formule. Integrac´ı kvadratick´eho interpolaˇcn´ıho polynomu P2(x) proch´azej´ıc´ıho body [a, f(a)], [1 2(a+b), f(1 2(a+b))] a [b, f(b)] dostaneme QS(f) = Zb a P2(x) dx=b−a 6f(a) + 4fa+b 2+f(b).(4.20) Simpsonova formule je ˇr´adu 3. Ovˇeˇrte! (Staˇc´ı porovnat QS(xk) a I(xk) na intervalu h0,1i postupnˇe pro k= 0,1,2,3,4.) Kdyˇz f∈C4ha, bi, pro chybu plat´ı RS(f) = −1 90f(4)(ξ)b−a 25 ,kde ξ∈(a, b).(4.21) Booleova formule vznikne integrac´ı interpolaˇcn´ıho polynomu P4(x), jehoˇz uzly jsou, kromˇe koncov´ych bod˚u a,b, tak´e stˇred 1 2(a+b) intervalu ha, bia body a+1 4(b−a) 70 ab−1 4(b−a) leˇz´ıc´ı v jedn´e ˇctvrtinˇe a ve tˇrech ˇctvrtin´ach intervalu ha, bi. Booleova formule QB(f) = b−a 90 7f(a) + 32f3a+b 4+ 12fa+b 2+ 32fa+ 3b 4+ 7f(b) je ˇr´adu 5 (ovˇeˇrte!). Pokud f∈C6ha, bi, pro chybu plat´ı RB(f) = −8 945f(6)(ξ)b−a 47 ,kde ξ∈(a, b). 4.3.2. Sloˇzen´e formule Abychom dostali dostateˇcnˇe pˇresnou aproximaci integr´alu I(f), rozdˇel´ıme interval ha, bina kratˇs´ı podintervaly a na kaˇzd´em z nich pouˇzijeme nˇekterou ze z´akladn´ıch formul´ı. Omez´ıme se na pˇr´ıpad, kdy z´akladn´ı formule na podintervalech jsou vˇzdy stejn´e, a to bud’to obd´eln´ıkov´e nebo lichobˇeˇzn´ıkov´e nebo Simpsonovy (sestaven´ı sloˇzen´e Booleovy formule ponech´av´ame ˇcten´aˇri jako cviˇcen´ı). Budeme uvaˇzovat rovnomˇern´e (= ekvidistantn´ı) dˇelen´ı a=x0< x1<···< xn=b, kde xi=a+ih,h= (b−a)/n ai= 0,1,...,n. (4.22) Sloˇzenou formuli na dˇelen´ı (4.22) budeme znaˇcit pomoc´ı horn´ıho indexu n. D´elka hse naz´yv´a krok dˇelen´ı. Sloˇzenou obd´eln´ıkovou formuli dostaneme souˇctem jednoduch´ych obd´eln´ıkov´ych formul´ı hf(xi−1+1 2h) na podintervalech hxi−1, xii. V´ysledkem je sloˇzen´a formule Qn M(f) := hfx0+1 2h+fx1+1 2h+···+fxn−1+1 2h.(4.23) Chybu Rn M(f) dostaneme jako souˇcet d´ılˇc´ıch chyb 1 24f′′(ξi)h3na hxi−1, xii: Rn M(f) = 1 24h3f′′(ξ1) + 1 24h3f′′(ξ2) + ···+1 24h3f′′(ξn) = 1 24h2b−a n[f′′(ξ1) + f′′(ξ2) + ···+f′′(ξn)] . Ze spojitosti f′′ plyne, ˇze aritmetick´y pr˚umˇer [f′′(ξ1) + f′′(ξ2) + ···+f′′(ξn)]/n druh´ych derivac´ı je roven druh´e derivacif′′(ξ) v nˇejak´em bodˇe ξ∈(a, b). Pro chybu Rn M(f) sloˇzen´e obd´eln´ıkov´e formule tedy plat´ı Rn M(f) = b−a 24 f′′(ξ)h2,kde ξ∈(a, b).(4.24) Sloˇzen´a lichobˇeˇzn´ıkov´a formule vznikne souˇctem jednoduch´ych lichobˇeˇzn´ıkov´ych formul´ı 1 2h[f(xi−1) + f(xi)] na podintervalech hxi−1, xii. V´ysledkem je formule Qn T(f) := h1 2f(x0) + f(x1) + ···+f(xn−1) + 1 2f(xn).(4.25) 71 x0x1x2xi−1 xixi+1 xn x0x1x2xi−1 xixi+1 xn Obr. 4.3: Sloˇzen´a obd´eln´ıkov´a a lichobˇeˇzn´ıkov´a formule Chybu Rn T(f) dostaneme jako souˇcet chyb −1 12 f′′(ξi)h3na podintervalech hxi−1, xii. Po jednoduch´e ´upravˇe obdrˇz´ıme Rn T(f) = −b−a 12 f′′(ξ)h2,kde ξ∈(a, b).(4.26) V MATLABu lze pouˇz´ıt funkci trapz. Sloˇzenou Simpsonovu formuli dostaneme pro sud´y poˇcet nd´ılk˚u tak, ˇze seˇcteme jednoduch´e Simpsonovy formule na intervalech hx0, x2i,hx2, x4i,...,hxn−2, xnid´elky 2h: Qn S(f) :=2h 6[f(x0) + 4f(x1) + f(x2)] + 2h 6[f(x2) + 4f(x3) + f(x4)] + ···+ 2h 6[f(xn−4) + 4f(xn−3) + f(xn−2)] + 2h 6[f(xn−2) + 4f(xn−1) + f(xn)] . Odtud tedy Qn S(f) := h 3[f(x0)+4f(x1)+2f(x2)+4f(x3)+···+2f(xn−2)+4f(xn−1)+f(xn)] .(4.27) 72 Chyba Rn S(f) je souˇctem chyb na podintervalech: Rn S(f) = −1 90h5f(4)(ξ2)−1 90h5f(4)(ξ4)−...−1 90h5f(4)(ξn) = −1 90h4b−a 22 nf(4)(ξ2) + f(4)(ξ4) + ···+f(4)(ξn). Kdyˇz aritmetick´y pr˚umˇer ˇctvrt´ych derivac´ı (tj. v´yraz ve sloˇzen´e z´avorce) nahrad´ıme ˇclenem f(4)(ξ), dostaneme Rn S(f) = −b−a 180 f(4)(ξ)h4,kde ξ∈(a, b).(4.28) Jak dos´ahnout poˇzadovan´e pˇresnosti. Vyvst´av´a ot´azka jak zajistit, aby chyba, kter´e se pˇri numerick´em v´ypoˇctu integr´alu dopust´ıme, byla menˇs´ı neˇz zadan´a tolerance ε. Ve vzorc´ıch (4.24), (4.26) a (4.28) bohuˇzel vystupuje bl´ıˇze neurˇcen´e ˇc´ıslo ξ, o nˇemˇz v´ıme jen to, ˇze leˇz´ı nˇekde v intervalu (a, b). Kdyˇz oznaˇc´ıme M2= max x∈ha,bi|f′′(x)|, M4= max x∈ha,bi|f(4)(x)|, pak chyby m˚uˇzeme odhadnout podle vztahu |Rn M(f)| ≤ b−a 24 M2h2pro sloˇzenou obd´eln´ıkovou formuli, (4.24’) |Rn T(f)| ≤ b−a 12 M2h2pro sloˇzenou lichobˇeˇzn´ıkovou formuli, (4.26’) |Rn S(f)| ≤ b−a 180 M4h4pro sloˇzenou Simpsonovu formuli. (4.28’) Poˇcet d´ılk˚u n= (b−a)/h pak lze urˇcit tak, aby |Rn M(f)| ≤ εpopˇr. |Rn T(f)| ≤ εnebo |Rn S(f)| ≤ ε. Takto stanoven´y poˇcet d´ılk˚u je vˇsak zpravidla pˇrehnanˇe velk´y. To je d˚usledek toho, ˇze jsme v odhadech chyb nahradili derivaci v bl´ıˇze neurˇcen´em bodˇe maximem t´eto derivace na cel´em intervalu ha, bi. Pˇr´ıklad 4.3. Urˇceme poˇcet d´ılk˚u ntak, abychom vypoˇc´ıtali Rπ/2 0excos xdxs chybou nejv´yˇse 10−4pomoc´ı sloˇzen´e obd´eln´ıkov´e, lichobˇeˇzn´ıkov´e a Simpsonovy formule. Pˇresn´a hodnota Rπ/2 0excos xdx=1 2ex(cos x+ sin x)π/2 0 . = 1,905239. Pro f(x) = excos xje f′′(x) = −2exsin xa tedy M2≤2eπ/2. V pˇr´ıpadˇe obd´eln´ıkov´e formule pomoc´ı (4.24’) odvod´ıme: |Rn M(f)| ≤ π/2 24 2eπ/2π/2 n2 ≤10−4=⇒n≥125 . Provedeme-li v´ypoˇcet s t´ımto npodle (4.23), dostaneme Q125 M(f). = 1,905277, takˇze skuteˇcn´a chyba je asi 3,8·10−5. Cestou pokus˚u se uk´azalo, ˇze pro dosaˇzen´ı poˇzadovan´e tolerance staˇc´ı vz´ıt n= 78. 73 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −6 −4 −2 0 2 4 x y y = 4sinx−x3−1 Obr. 5.1: Graf funkce 4 sin x−x3−1 Stˇred xk+1 intervalu (ak, bk) aproximuje koˇren x∗s chybou |xk+1 −x∗| ≤ 1 2(bk−ak) = 2−k−1(b0−a0).(5.2) Pro k→ ∞ zˇrejmˇe |Ik| → 0 a xk→x∗. Pˇr´ıklad 5.2. Metodu bisekce aplikujeme na rovnici z pˇr´ıkladu 5.1. Jako poˇc´ateˇcn´ı zvok akbkxk+1 f(xk+1) 0 1 2 1,5<0 1 1 1,5 1,25 >0 2 1,25 1,5 1,375 >0 3 1,375 1,5 1,4375 <0 4 1,375 1,4375 1,40625 >0 5 1,40625 1,4375 1,421875 l´ıme interval (a0, b0) = (1,2). Pˇripomeˇnme, ˇze f(1) >0, f(2) <0. Proto tak´e f(ak)≥0, f(bk)≤0 pro kaˇzd´e k. Posloupnost interval˚u zaznamen´av´ame do tabulky. Po pˇeti kroc´ıch m´a interval (a5, b5) d´elku 2−5= 0,03125 a x6= 1,421875 aproximuje koˇren s chybou nepˇresahuj´ıc´ı 2−6= 0,015625.  Metoda bisekce konverguje pomalu: protoˇze 10−1. = 2−3,32, zpˇresnˇen´ı o jednu dekadickou cifru vyˇzaduje v pr˚umˇeru 3,32 kroku. Vˇsimnˇete si, ˇze rychlost konvergence vyj´adˇren´a vztahem (5.2) v˚ubec nez´avis´ı na funkci f(x). To proto, ˇze jsme vyuˇz´ıvali pouze znam´enka funkˇcn´ıch hodnot. Kdyˇz tyto hodnoty (a pˇr´ıpadnˇe tak´e hodnoty derivac´ı f(x)) vyuˇzijeme efektivnˇeji, m˚uˇzeme dos´ahnout podstatnˇe rychlejˇs´ı konvergence. Takov´e zpˇresˇnuj´ıc´ı metody vˇsak konverguj´ı pouze tehdy, kdyˇz pro nˇe zvol´ıme dostateˇcnˇe dobrou poˇc´ateˇcn´ı aproximaci. Vhodn´a poˇc´ateˇcn´ı aproximace b´yv´a ˇcasto urˇcena pr´avˇe metodou bisekce. 5.2. Zpˇresˇnuj´ıc´ı metody Snad nejzn´amˇejˇs´ı mezi nimi je Newtonova metoda nebo-li metoda teˇcen. Jak je u iteraˇcn´ıch metod zvykem, vyjdeme z poˇc´ateˇcn´ı aproximace x0a postupnˇe poˇc´ıt´ame x1, x2,... zp˚usobem, kter´y si ted’vysvˇetl´ıme. Pˇredpokl´adejme, ˇze zn´ame xka m´ame urˇcit lepˇs´ı aproximaci xk+1. Udˇel´ame to tak, ˇze bodem [xk, f(xk)] vedeme teˇcnu ke kˇrivce y=f(x) a pr˚useˇc´ık teˇcny s osou xpovaˇzujeme 80 za xk+1. Do rovnice teˇcny y=f(xk) + f′(xk)(x−xk) tedy dosad´ıme y:= 0, vypoˇcteme xa poloˇz´ıme xk+1 := x. Tak dostaneme pˇredpis xk+1 =xk−f(xk) f′(xk).(5.3) V´ypoˇcet ukonˇc´ıme a xk+1 povaˇzujeme za dostateˇcnˇe pˇresnou aproximaci koˇrene, pokud |xk+1 −xk| ≤ ε , pˇr´ıpadnˇe |xk+1 −xk| ≤ ε|xk|nebo |f(xk+1)| ≤ ε , (5.4) kde εje poˇzadovan´a pˇresnost. T´ım sice nen´ı zaruˇceno, ˇze tak´e |xk+1 −x∗| ≤ ε, je to ale obvykl´y zp˚usob, pomoc´ı nˇehoˇz iterace ukonˇc´ıme. Tato tzv. stop krit´eria jsou vhodn´a i pro dalˇs´ı metody, kter´e v tomto odstavci uvedeme. x0 x1 x2 x3 f(x) y=0 [xi,f(xi)] (y−yi)/(x−xi) = f ’(xi) Obr. 5.2: Newtonova metoda Pˇr´ıklad 5.3. Newtonovou metodou urˇc´ıme kladn´y koˇren rovnice z pˇr´ıkladu 5.1. Zvol´ıme x0= 2. V´ypoˇcet ukonˇc´ıme, kdyˇz |f(xk)|<10−5. Posledn´ı sloupec vyˇzaduje znalost pˇresk xkf(xk)f′(xk)f(xk)/f′(xk)xk−x∗ 0 2 −5,362810 −13,66459 0,392460 0,563550 1 1,607540 −1,156877 −7,899490 0,146450 0,171089 2 1,461090 −0,143158 −5,966406 0,023994 0,024640 3 1,437096 −0,003653 −5,662524 0,000645 0,000646 4 1,436451 −0,000003 0,000000 n´eho ˇreˇsen´ı. To z´ısk´ame proveden´ım jeˇstˇe jednoho kroku Newtonovy metody. D´a se uk´azat, ˇze x∗. =x5= 1,43645032 m´a vˇsechny cifry platn´e. Poˇzadovan´a pˇresnost byla tedy dosaˇzena ve ˇctvrt´em kroku, x4. = 1,43645 m´a vˇsechny cifry platn´e.  81 Konvergence Newtonovy metody. Necht ’ek=xk−x∗je chyba v k-t´em kroku. Uk´aˇzeme si, jak souvis´ı s chybou ek+1 v kroku n´asleduj´ıc´ım. Z Taylorova rozvoje f(x∗) okolo xkdostaneme 0 = f(x∗) = f(xk) + (x∗−xk)f′(xk) + 1 2(x∗−xk)2f′′(ξ), kde ξje nˇejak´y bl´ıˇze neurˇcen´y bod intervalu, jehoˇz krajn´ı body jsou xkax∗. Kdyˇz rovnici dˇel´ıme f′(xk), dostaneme −1 2(x∗−xk)2f′′(ξ) f′(xk)=f(xk) f′(xk)+ (x∗−xk) = x∗−xk−f(xk) f′(xk)=x∗−xk+1 , takˇze m´ame ek+1 =1 2 f′′(ξ) f′(xk)e2 k,(5.5) a kdyˇz xk→x∗, pak ek+1 e2 k−→ C , kde C=1 2 f′′(x∗) f′(x∗). Protoˇze chyba ek+1 je ´umˇern´a druh´e mocninˇe chyby ek, ˇr´ık´ame, ˇze Newtonova metoda konverguje kvadraticky nebo tak´e, ˇze je druh´eho ˇr´adu. Uved’me si pˇresnˇejˇs´ı definici: Necht ’x0, x1, x2,... je posloupnost, kter´a konverguje k x∗aek=xk−x∗. Kdyˇz existuje ˇc´ıslo pa konstanta C6= 0 takov´a, ˇze lim k→∞ |ek+1| |ek|p=C , (5.6) pak pse naz´yv´a ˇr´ad konvergence posloupnosti a Cje chybov´a konstanta. Speci´alnˇe ˇr´ık´ame, ˇze konvergence je line´arn´ı, superline´arn´ı, kvadratick´a, kdyˇz p= 1 aC < 1, p > 1, p= 2 . ˇ Rekneme, ˇze dan´a metoda je ˇr´adu p, jestliˇze vˇsechny konvergentn´ı posloupnosti z´ıskan´e touto metodou maj´ı ˇr´ad konvergence vˇetˇs´ı nebo rovn´y pa nejm´enˇe jedna z tˇechto posloupnost´ı m´a ˇr´ad konvergence rovn´y pˇresnˇe p. V bl´ızkosti koˇrene plat´ı: ˇc´ım vyˇsˇs´ı ˇr´ad p, t´ım rychlejˇs´ı konvergence, nebot ’ |ek+1| ≈ C|ek|p, takˇze kdyˇz |ek|je mal´e, pak |ek+1|je t´ım menˇs´ı, ˇc´ım je pvˇetˇs´ı. V´ıme uˇz, ˇze kdyˇz Newtonova metoda konverguje, pak rychlost konvergence xk→x∗ je alespoˇn kvadratick´a (pro nˇekter´e funkce fm˚uˇze b´yt i vyˇsˇs´ı). Zb´yv´a jeˇstˇe zodpovˇedˇet ot´azku, za jak´ych podm´ınek je zaruˇceno, ˇze konvergence v˚ubec nastane. Ukaˇzme si to. 82 Pˇredpokl´adejme, ˇze v nˇejak´em okol´ı Ikoˇrene plat´ı 1 2 f′′(y) f′(x)≤mpro vˇsechna x, y ∈I . Kdyˇz xk∈I, pak z (5.5) plyne |ek+1| ≤ m|ek|2nebo-li |mek+1| ≤ |mek|2. Opakov´an´ım t´eto ´uvahy dostaneme |mek+1| ≤ |mek|2≤ |mek−1|4≤ |mek−2|8≤ |mek−3|16 ≤ ··· ≤ |me0|r,kde r= 2k+1. Kdyˇz plat´ı |me0|<1, pak jistˇe |ek+1| → 0 a tedy xk+1 →x∗. Dok´azali jsme tedy, ˇze Newtonova metoda vˇzdy konverguje za pˇredpokladu, ˇze poˇc´ateˇcn´ı aproximaci zvol´ıme dostateˇcnˇe bl´ızko ke koˇrenu. Dobrou poˇc´ateˇcn´ı aproximaci x0m˚uˇzeme z´ıskat napˇr. metodou bisekce. Vhodn´ym spojen´ım metody bisekce a Newtonovy metody lze sestrojit kombinovanou metodu, kter´a vˇzdy konverguje, viz napˇr. procedura rtsafe v [18]. V bl´ızkosti koˇrene se pˇritom uplatn´ı jen Newtonova metoda, takˇze konvergence je rychl´a. Pomoc´ı n´aˇcrtku snadno ovˇeˇr´ıme, ˇze Newtonova metoda konverguje, kdyˇz jsou splnˇeny tzv. Fourierovy podm´ınky: a) f∈C2ha, bia pˇritom f(a)f(b)<0; b) f′af′′ nemˇen´ı na intervalu ha, biznam´enko a f′(x)6= 0 pro kaˇzd´e x∈ ha, bi; c) jako x0vol´ıme ten z bod˚u a, b, v nˇemˇz je f(x0)f′′(x0)>0. Praktick´y v´yznam vˇsak Fourierovy podm´ınky nemaj´ı, nebot ’pro velk´e b−aobvykle tyto podm´ınky bud’to neplat´ı nebo je neum´ıme snadno ovˇeˇrit. Metoda seˇcen. V kaˇzd´em kroku Newtonovy metody mus´ıme poˇc´ıtat hodnotu f(xk) a derivaci f′(xk). Kdyˇz vzorec pro v´ypoˇcet derivace nem´ame k dispozici, nebo kdyˇz n´aklady spojen´e s v´ypoˇctem derivace jsou vysok´e, m˚uˇzeme derivaci aproximovat pod´ılem f′(xk)≈f(xk)−f(xk−1) xk−xk−1 . Tak dostaneme metodu seˇcen: zad´ame dvˇe poˇc´ateˇcn´ı aproximace x0,x1a poˇc´ıt´ame x2, x3,... podle pˇredpisu xk+1 =xk−xk−xk−1 f(xk)−f(xk−1)f(xk).(5.7) N´azev metody vych´az´ı z jej´ı geometrick´e interpretace: xk+1 je x-ov´a souˇradnice pr˚useˇc´ıku pˇr´ımky proch´azej´ıc´ı body [xk−1, f(xk−1)] a [xk, f(xk)] s osou x: y=f(xk) + f(xk)−f(xk−1) xk−xk−1 (x−xk) = 0 =⇒x=xk+1 . Protoˇze tato pˇr´ımka prot´ın´a graf funkce f, je to seˇcna, odtud metoda seˇcen. 83 x0 x1 x2 x3 x4 f(x) y=0 [xi,f(xi)] (y−yi)/(x−xi) = (yi−yi−1)/(xi−xi−1) Obr. 5.3: Metoda seˇcen Vˇsimnˇete si, ˇze v kaˇzd´em kroku vyˇc´ıslujeme hodnotu funkce jen jednou: vypoˇcteme f(xk), hodnotu f(xk−1) pˇrevezmeme z pˇredchoz´ıho kroku. D´a se odvodit, ˇze rychlost konvergence metody seˇcen je ˇr´adu p=1 2(1 + √5) ≈1,618, tedy ponˇekud niˇzˇs´ı neˇz u Newtonovy metody. ˇ C´ıslo τ= (√5−1)/2≈0,618 je tzv. pomˇer zlat´eho ˇrezu. Toto magick´e ˇc´ıslo se znovu objev´ı v odstavci 6.1, kde si o nˇem ˇrekneme trochu v´ıc (viz pozn´amka o zlat´em ˇrezu). Pˇr´ıklad 5.4. Metodou seˇcen urˇc´ıme kladn´y koˇren rovnice z pˇr´ıkladu 5.1. Zvol´ıme x0= 1, k xkf(xk)xk−x∗ 0 1 1,365884 −0,436450 1 2 −5,362810 0,563550 2 1,202994 0,991513 −0,233456 3 1,327357 0,543420 −0,109094 4 1,478177 −0,246970 0,041726 5 1,431051 0,030349 −0,005400 6 1,436208 0,001370 −0,000242 7 1,436452 −0,000008 0,000001 x1= 2. V´ypoˇcet ukonˇc´ıme, kdyˇz bude |f(xk)|<10−5. Aˇz do ˇctvrt´eho kroku (v´ypoˇcet x5) je konvergence pomˇernˇe pomal´a. Teprve v posledn´ıch dvou kroc´ıch se plnˇe uplatnila rychl´a konvergence metody seˇcen.  Metoda seˇcen zaruˇcenˇe konverguje, pokud zvol´ıme startovac´ı hodnoty x0ax1dostateˇcnˇe bl´ızko ke koˇrenu x∗. To lze zajistit napˇr. metodou bisekce. Dalˇs´ı metodou, jak z´ıskat dobr´e startovac´ı aproximace, je varianta metody seˇcen zn´am´a jako Metoda regula falsi. Poˇc´ateˇcn´ı aproximace x0ax1se vol´ı tak, aby f(x0)f(x1)<0. Nov´a aproximace xk+1 se opˇet z´ısk´a jako pr˚useˇc´ık seˇcny s osou x. Seˇcna vˇsak tentokr´at spojuje bod [xk, f(xk)] s bodem [xℓ, f(xℓ)], kde ℓje nejvˇetˇs´ı index, pro kter´y f(xk)f(xℓ)<0. V´ypoˇcet tedy prob´ıh´a podle vzorce xk+1 =xk−xk−xℓ f(xk)−f(xℓ)f(xk), k = 1,2,.... (5.8) 84 Pˇritom pro k= 1 je ℓ= 0, a po v´ypoˇctu xk+1 urˇc´ıme index ℓtakto: kdyˇz f(xk+1)f(xℓ)>0, pak ℓ=k, v opaˇcn´em pˇr´ıpadˇe se ℓnemˇen´ı. V´yhodou metody regula falsi je to, ˇze podobnˇe jako metoda bisekce vˇzdy konverguje: interval Ik, jehoˇz koncov´e body jsou xkaxℓ, obsahuje koˇren. Na rozd´ıl od metody bisekce vˇsak d´elka intervalu Iknekonverguje k nule. Rychlost konvergence metody regula falsi je jen line´arn´ı. Metodu regula falsi (podobnˇe jako metodu bisekce) proto pouˇz´ıv´ame pouze pro z´ısk´an´ı dobr´e poˇc´ateˇcn´ı aproximace, pak pˇrech´az´ıme na rychlejˇs´ı metodu. x0 x1 x2 x3 x4 Obr. 5.4: Regula falsi Pˇr´ıklad 5.5. Metodou regula falsi urˇc´ıme kladn´y koˇren rovnice z pˇr´ıkladu 5.1. Zvol´ıme k ℓ xℓxkf(xk)xk−x∗ 0 1 1,365884 −0,436450 1 0 1 2 −5,362810 0,563550 2 1 2 1,202994 0,991513 −0,233456 3 1 2 1,327357 0,543420 −0,109094 4 1 2 1,389245 0,253012 −0,047205 5 1 2 1,416762 0,108896 −0,019688 6 1 2 1,428369 0,045283 −0,008081 7 1 2 1,433156 0,018561 −0,003295 . . . 15 1 2 1,436448 0,000014 −0,000002 16 1 2 1,436449 0,000006 −0,000001 x0= 1, x1= 2. Z tabulky je vidˇet, ˇze poˇc´ınaje druh´ym krokem je xℓ= 2. D´ale je zˇrejm´e, ˇze do ˇctvrt´eho kroku je pˇresnost metody regula falsi srovnateln´a s pˇresnost´ı metody seˇcen, viz pˇr´ıklad 5.4. V n´asleduj´ıc´ıch kroc´ıch je uˇz ale patrn´a line´arn´ı konvergence, podm´ınka |f(xk)|<10−5je splnˇena aˇz pro x16. Vˇsimnˇete si: d´elka interval˚u Ik= (xk,2), k≥2, konverguje k ˇc´ıslu x−x∗. = 0,563550.  Steffensenova metoda se ˇr´ıd´ı pˇredpisem xk+1 =xk−f(xk) dk ,kde dk=f(xk+f(xk)) −f(xk) f(xk)(5.9) 85 je speci´alnˇe spoˇcten´a aproximace f′(xk) pˇripom´ınaj´ıc´ı dopˇrednou diferenci: f′(xk)≈dk=f(xk+hk)−f(xk) hk ,kde hk=f(xk). V kaˇzd´em kroku se funkce fvyhodnocuje dvakr´at: kromˇe hk=f(xk) se poˇc´ıt´a jeˇstˇe tak´e f(xk+hk). Oproti metodˇe seˇcen je tu jedno vyhodnocen´ı funkce nav´ıc. Na druh´e stranˇe lze uk´azat, ˇze rychlost konvergence Steffensenovy metody je stejn´a jako u Newtonovy metody, tedy kvadratick´a. Metoda inverzn´ı kvadratick´e interpolace. Metoda seˇcen pouˇz´ıv´a dva pˇredchoz´ı body k z´ısk´an´ı dalˇs´ıho, proˇc tedy nepouˇz´ıt tˇri? Body [xk−2, f(xk−2)], [xk−1, f(xk−1)] a [xk, f(xk)] m˚uˇzeme proloˇzit parabolu P2(x) a hledat jej´ı pr˚useˇc´ık s osou x. Za dalˇs´ı aproximaci xk+1 pak zvol´ıme ten z koˇren˚u polynomu P2(x), kter´y je bl´ıˇz k pˇredchoz´ı aproximaci xk. Na tomto principu je zaloˇzena M˝ullerova metoda. Pot´ıˇz je v tom, ˇze parabola nemus´ı x-ovou osu protnout, nebot ’kvadratick´a funkce P2(x) nemus´ı m´ıt re´aln´e koˇreny. V´ypoˇcet je proto tˇreba prov´adˇet v komplexn´ı aritmetice, a to i v pˇr´ıpadˇe, ˇze rovnice f(x) = 0 m´a jen re´aln´e koˇreny. M´ısto paraboly v promˇenn´e xm˚uˇzeme tˇremi body proloˇzit parabolu Q2(y) v promˇenn´e y, urˇcenou interpolaˇcn´ımi podm´ınkami Q2(f(xk−2)) = xk−2, Q2(f(xk−1)) = xk−1, Q2(f(xk)) = xk. Jsou-li hodnoty f(xk−2), f(xk−1) a f(xk) navz´ajem r˚uzn´e, parabola Q2(y) existuje a prot´ın´a osu xv jedin´em bodˇe. Klademe tedy xk+1 =Q2(0). Tato metoda je zn´ama jako metoda inverzn´ı kvadratick´e interpolace. Jej´ı konvergence je superline´arn´ı ˇr´adu p≈1,839, viz [7]. Brentova metoda. Metoda inverzn´ı kvadratick´e interpolace spolu s metodou seˇcen a metodou bisekce jsou z´akladem popul´arn´ı Brentovy metody, viz napˇr. [15], d´ale tak´e program zbrent v [18] nebo funkce fzero v MATLABu. Pˇrednost´ı Brentovy metody je to, ˇze nepouˇz´ıv´a derivace funkce f, je spolehliv´a, tj. zaruˇcenˇe konverguje ke koˇrenu, a po nˇekolika poˇc´ateˇcn´ıch kroc´ıch se chyba rychle zmenˇsuje, nebot ’rychlost konvergence je superline´arn´ı. Startovac´ı body x0ax1je tˇreba zvolit tak, aby f(x0)f(x1)<0. Aproximace x2se urˇc´ı metodou seˇcen. Necht ’(a1, b1) je interval, jehoˇz koncov´e body jsou x0ax1. Pak zˇrejmˇe x2∈(a1, b1). Dalˇs´ı aproximaci x3budeme hledat v kratˇs´ım intervalu (a2, b2)⊂(a1, b1), jehoˇz jeden koncov´y bod je x2a druh´y je ten z bod˚u a1,b1, v nˇemˇz m´a funkce fopaˇcn´e znam´enko neˇz v x2, takˇze f(a2)f(b2)<0a(a2, b2) obsahuje koˇren. Pˇri v´ypoˇctu x3, x4,... Brentova metoda pouˇz´ıv´a jednu ze tˇr´ı z´akladn´ıch metod tak, aby nov´a aproximace xk+1 ∈(ak, bk). D´ale se vybere interval (ak+1, bk+1)⊂(ak, bk) obsahuj´ıc´ı koˇren. Jedn´ım z jeho koncov´ych bod˚u je xk+1, druh´ym je ten z bod˚u ak,bk, v nˇemˇz m´a funkce fznam´enko opaˇcn´e neˇz v xk+1. Pˇri v´ypoˇctu xk+1 se pˇrednostnˇe pouˇzije metoda inverzn´ı kvadratick´e interpolace, pokud takto z´ıskan´a aproximace nen´ı dostateˇcnˇe dobr´a, zkus´ı se metoda seˇcen, a kdyˇz ani ta nezabere, pouˇzije se jako z´achrana metoda bisekce. Podrobnˇejˇs´ı popis Brentovy metody je uveden napˇr. v [15], [18]. Pˇr´ıklad 5.6. Budeme hledat kladn´y koˇren rovnice z pˇr´ıkladu 5.1 a porovn´ame jednotliv´e metody podle poˇctu pk krok˚u a poˇctu pf vyhodnocen´ı funkce f(u Newtonovy metody 86 do pf zahrneme tak´e poˇcet vyhodnocen´ı derivace f′). Pro v´ypoˇcet Brentovou metodou jsme pouˇzili upraven´y program fzerotx, viz [15]. V´ypoˇcet jsme zah´ajili takto: v metodˇe bisekce poˇc´ateˇcn´ı interval (a0, b0) = (1,2), v Newtonovˇe a Steffensenovˇe metodˇe x0= 2, v ostatn´ıch metod´ach x0= 1 a x1= 2. Pouˇzili jsme stop krit´erium |f(xk)|< ε. V tabulce jsou uvedeny hodnoty pk/pf pro nˇekolik toleranc´ı ε. ε10−310−610−910−12 10−15 bisekce 9/11 19/21 29/31 39/41 49/51 regula falsi 10/12 17/19 25/27 33/35 40/42 seˇcny 6/8 7/9 8/10 8/10 9/11 Newton 4/8 5/10 5/10 6/12 6/12 Steffensen 4/8 5/10 6/12 6/12 7/14 Brent 6/7 7/8 7/8 8/9 8/9 Nejmenˇs´ı pk m´a Newtonova metoda, nejmenˇs´ı pf Brentova metoda. Z v´ypisu o pr˚ubˇehu v´ypoˇctu Brentovou metodou vypl´yv´a, ˇze se ani jednou nepouˇzila bisekce, proto tak skvˇel´y v´ysledek.  Pozn´amka (O metodˇe prost´e iterace). Pˇredpokl´adejme, ˇze funkce g∈Cha, bisplˇnuje tyto dvˇe podm´ınky: (α)g(x)∈ ha, bi ∀x∈ ha, bi, (β) existuje ˇc´ıslo q, 0 ≤q < 1, takov´e, ˇze |g(x)−g(y)| ≤ q|x−y| ∀x, y ∈ ha, bi. Pak rovnice x=g(x) m´a v ha, bijedin´e ˇreˇsen´ı x∗aposloupnost postupn´ych aproximac´ı xk+1 =g(xk), k= 0,1,..., konverguje k x∗pro kaˇzd´e x0∈ ha, bi. Bod x∗=g(x∗) se naz´yv´a pevn´y bod funkce g(zobrazuje x∗na sebe). N´asleduje n´aˇcrt d˚ukazu. 1) Existence. Z podm´ınky (α) plyne g(a)≥a,g(b)≤b, odtud (a−g(a))·(b−g(b)) ≤0, takˇze v ha, bileˇz´ı koˇren rovnice x−g(x) = 0. 2) Jednoznaˇcnost. Necht ’pro x∗, y∗∈ ha, biplat´ı x∗=g(x∗), y∗=g(y∗). Podle (β) |x∗−y∗|=|g(x∗)−g(y∗)| ≤ q|x∗−y∗|, coˇz je moˇzn´e jedinˇe kdyˇz x∗=y∗. 3) Konvergence. Podle (β) je |xk−x∗|=|g(xk−1)−g(x∗)| ≤ q|xk−1−x∗|. Opakov´an´ım t´eto ´uvahy dostaneme nakonec |xk−x∗| ≤ qk|x0−x∗| → 0 pro k→ ∞, takˇze xk→x∗. M´ısto podm´ınky (β) m˚uˇzeme pro g∈C1ha, bipouˇz´ıt silnˇejˇs´ı podm´ınku (β′)|g′(x)| ≤ q < 1∀x∈ ha, bi. Podle vˇety o stˇredn´ı hodnotˇe totiˇz g(x)−g(y) = g′(ξ)(x−y), kde ξleˇz´ı mezi xay, takˇze pro x, y ∈ ha, bipodle (β′) je |g(x)−g(y)|=|g′(ξ)|·|x−y| ≤ q|x−y|, tj. plat´ı (β). Vˇsimnˇete si, ˇze pro ha, bi=hx∗−δ, x∗+δije platnost podm´ınky (α) d˚usledkem platnosti podm´ınky (β): |g(x)−x∗|=|g(x)−g(x∗)| ≤ q|x−x∗|<|x−x∗|, tj. kdyˇz |x−x∗| ≤ δ, pak tak´e |g(x)−x∗| ≤ δ. Pˇribliˇzn´y v´ypoˇcet koˇrene x∗rovnice x=g(x) podle formule xk+1 =g(xk) se naz´yv´a metoda prost´e iterace nebo tak´e metoda postupn´ych aproximac´ı. Podm´ınky (α) a (β) nebo (β′) jsou postaˇcuj´ıc´ı pro konvergenci t´eto metody. 87 Vhodnou ´upravou rovnice f(x) = 0 na tvar x=g(x) m˚uˇzeme dostat ˇradu r˚uzn´ych konkr´etn´ıch metod. Tak tˇreba pro g(x) = x−f(x)/f′(x) dostaneme Newtonovu metodu. Rychlost konvergence posloupnosti postupn´ych aproximac´ı {xk}∞ k=0 z´avis´ı na chov´an´ı funkce gv bodˇe x∗. Jsou-li splnˇeny podm´ınky (α) a (β) nebo (β′), a m´a-li gdostateˇcn´y poˇcet spojit´ych derivac´ı, daj´ı se dok´azat n´asleduj´ıc´ı tvrzen´ı. •Pokud g′(x∗)6= 0, je ˇr´ad konvergence roven jedn´e a plat´ı |xk+1 −x∗| ≤ q|xk−x∗|. •Pokud g′(x∗) = 0 a g′′(x∗)6= 0, je ˇr´ad konvergence roven dvˇema. •Obecnˇe, pokud jsou derivace g(s)(x∗) = 0, s= 1,2,...,r−1, a g(r)(x∗)6= 0, konvergence je ˇr´adu r. Pro Newtonovu metodu g′(x) = f(x)f′′(x)/(f′(x))2, tj. g′(x∗) = 0, takˇze konvergence xk→x∗je ˇr´adu alespoˇn dva (coˇz potvrzuje n´am jiˇz zn´am´y v´ysledek). D´a se tak´e dok´azat, ˇze kdyˇz |g′(x∗)|>1, pak pro x06=x∗posloupnost postupn´ych aproximac´ı k x∗konvergovat nem˚uˇze.  Pˇr´ıklad 5.7. Neline´arn´ı rovnice f(x) = x2−2x−3 = 0 m´a koˇreny x∗=−1 a x∗= 3. Prozkoum´ame konvergenci ke koˇrenu x∗= 3 pro nˇekolik iteraˇcn´ıch funkc´ı g. 1. g(x) = (x2−3)/2, g′(x) = x,|g′(3)|= 3, pro x06= 3 konvergence nenastane. 2. g(x) = √2x+ 3, g′(x) = 1/√2x+ 3, |g′(3)|= 1/3, line´arn´ı konvergence nastane napˇr. pro libovoln´e x0z intervalu h2,4i, nebot ’v nˇem |g′(x)| ≤ 1/√7. 3. g(x) = 2 + 3/x,g′(x) = −3/x2,|g′(3)|= 1/3, line´arn´ı konvergence nastane napˇr. pro libovoln´e x0z intervalu h2,4i, nebot ’v nˇem |g′(x)| ≤ 3/4. 4. g(x) = (x2+3)/(2x−2), g′(x) = 2(x2−2x−3)/(2x−2)2,g′(3) = 0, g′′(3) = 1/2, kvadratick´a konvergence nastane napˇr. pro x0z intervalu h2,5; 3,5i, v nˇemˇz |g′(x)|<0,39, jak snadno zjist´ıme. Ovˇeˇrte, ˇze tato iteraˇcn´ı funkce odpov´ıd´a Newtonovˇe metodˇe.  Pozn´amka (O n´asobn´ych koˇrenech). ˇ Rekneme, ˇze koˇren x∗rovnice f(x) = 0 m´a n´asobnost q, jestliˇze funkce g(x) = f(x)/(x−x∗)qje v bodˇe x∗definov´ana a koˇren v nˇem uˇz nem´a, tj. kdyˇz 0 <|g(x∗)|<∞. Jestliˇze m´a funkce f(x) v okol´ı koˇrene x∗spojit´e derivace aˇz do ˇr´adu qvˇcetnˇe, pak f(j)(x∗) = 0, j= 0,1,...,q−1. Nˇekter´e z doposud uveden´ych metod lze pouˇz´ıt tak´e pro nalezen´ı n´asobn´ych koˇren˚u, konvergence vˇsak b´yv´a pomalejˇs´ı. Tak tˇreba Newtonova metoda konverguje jen line´arnˇe s chybovou konstantou C= (q−1)/q. Kdyˇz oˇcek´av´ame, ˇze rovnice f(x) = 0 m˚uˇze m´ıt n´asobn´e koˇreny, je vhodn´e vyuˇz´ıt toho, ˇze funkce u(x) = f(x)/f′(x) m´a pouze jednoduch´e koˇreny. M´ısto rovnice f(x) = 0 tedy ˇreˇs´ıme rovnici u(x) = 0.  Pozn´amka (O dosaˇziteln´e pˇresnosti). Necht ’xkje aproximace jednoduch´eho koˇrene rovnice f(x) = 0. Pomoc´ı vˇety o stˇredn´ı hodnotˇe dostaneme f(xk) = f(xk)−f(x∗) = f′(ξ)(xk−x∗), kde ξje nˇejak´y bod leˇz´ıc´ı mezi xkax∗. Pˇredpokl´adejme, ˇze pˇri v´ypoˇctech pracujeme jen s pˇribliˇzn´ymi hodnotami ˜ f(xk) = f(xk) + δk, pˇriˇcemˇz |δk| ≤ δ. Pak nejlepˇs´ı v´ysledek, 88 kter´eho m˚uˇzeme dos´ahnout, je ˜ f(xk) = 0. V tom pˇr´ıpadˇe |f(xk)| ≤ δ, takˇze |xk−x∗|=|f(xk)| |f′(ξ)|≤δ |f′(ξ)|≈δ |f′(x∗)|=: ε∗ x, pokud se f′v bl´ızkosti koˇrene pˇr´ıliˇs nemˇen´ı. Vypoˇc´ıtat x∗s menˇs´ı chybou neˇz ε∗ xnelze. Proto se ε∗ xnaz´yv´a dosaˇziteln´a pˇresnost koˇrene x∗. Vˇsimnˇete si: kdyˇz je velikost smˇernice |f′(x∗)|v koˇrenu x∗mal´a, je dosaˇziteln´a pˇresnost ε∗ xvelk´a, viz obr. 5.5. V takov´em pˇr´ıpadˇe je v´ypoˇcet koˇrene x∗ˇspatnˇe podm´ınˇen´y probl´em: mal´a zmˇena fvyvol´a velkou zmˇenu x∗. x*−ε* xx*+ε* x f(x) y=0 Obr. 5.5: Dosaˇziteln´a pˇresnost koˇrene Podobn´a ´uvaha pro koˇren n´asobnosti qd´av´a dosaˇzitelnou pˇresnost ε∗ x=δ·q! f(q)(x∗)1/q . Exponent 1/q je pˇr´ıˇcinou toho, ˇze v´ypoˇcet n´asobn´eho koˇrene je obecnˇe ˇspatnˇe podm´ınˇen´a ´uloha. Tak tˇreba pro f(x) = xqje x∗= 0 koˇren n´asobnosti qaε∗ x=δ1/q. Pro q= 15 aδ= 10−15 dostaneme ε∗ x= 0,1!  Pozn´amka (O koˇrenech polynom˚u). Polynom pn(x) stupnˇe nm´a nobecnˇe komlexn´ıch koˇren˚u. Pro v´ypoˇcet jednoduch´ych re´aln´ych koˇren˚u funkce f(x) = pn(x) lze pouˇz´ıt libovolnou z dosud uveden´ych metod. O tom, jak se vypoˇr´adat s pˇr´ıpadn´ymi n´asobn´ymi koˇreny, pojedn´av´a v´yˇse uveden´a pozn´amka. Pro v´ypoˇcet komplexn´ıch koˇren˚u lze pouˇz´ıt napˇr. Newtonovu metodu, v n´ıˇz jako poˇc´ateˇcn´ı aproximaci vol´ıme komplexn´ı ˇc´ıslo. Pokud n´as zaj´ımaj´ı vˇsechny koˇreny polynomu, tak po nalezen´ı re´aln´eho koˇrene x∗ polynom pn(x) dˇel´ıme ˇclenem x−x∗. Tak dostaneme polynom pn−1(x) = pn(x)/(x−x∗) stupnˇe n−1 a d´ale hled´ame jeho koˇreny. Kdyˇz je x∗komplexn´ı koˇren, pak je koˇrenem tak´e komplexnˇe sdruˇzen´e ˇc´ıslo ¯x∗. V tom pˇr´ıpadˇe dˇel´ıme pn(x) kvadratick´ym polynomem (x−x∗)(x−¯x∗), jehoˇz koeficienty jsou re´aln´a ˇc´ısla. Tak dostaneme polynom pn−2(x) stupnˇe n−2 s re´aln´ymi koeficienty a pokraˇcujeme hled´an´ım jeho koˇren˚u. Pro v´ypoˇcet koˇren˚u polynom˚u jsou navrˇzeny tak´e speci´aln´ı, velmi efektivn´ı metody, o nichˇz lze z´ıskat informace napˇr. v [22].  5.3. Soustavy neline´arn´ıch rovnic Mnoh´e z metod urˇcen´ych pro ˇreˇsen´ı jedn´e neline´arn´ı rovnice lze zobecnit na ˇreˇsen´ı soustav neline´arn´ıch rovnic. Bohuˇzel to neplat´ı pro metodu bisekce ani pro metodu regula 89 5.4. Cviˇcen´ı 5.1. K jak´emu ˇc´ıslu konverguje posloupnost iterac´ı, definovan´a takto x0= 1, xk+1 =1 2xk+1 xk . Zd˚uvodnˇete konvergenci! [ lim k→∞ xk=√2,iteraˇcn´ı pˇredpis je Newtonova metoda pro rovnici x2−2 = 0. ] 5.2. Metodou seˇcen a teˇcen spoˇctˇete ˇreˇsen´ı d´ale uveden´ych rovnic na 6 platn´ych desetinn´ych cifer pˇresnˇe. Porovnejte poˇcet iterac´ı nutn´ych k dosaˇzen´ı poˇzadovan´e pˇresnosti. (a) x−e−x= 0, x0= 0,5, x1= 0,6; (b) x−cos x= 0, x0= 0,5, x1= 0,6; (c) x3+ 4x2−10 = 0, x0= 1,5, x1= 1,6. (x1pouˇzijte jen v metodˇe seˇcen.) [ (a) 0,567143 (b) 0,739085 (c) 1,365230.] 5.3. (a) Naprogramujte hled´an´ı znam´enkov´e zmˇeny. Vstupem je funkce f(x), interval ha, bia pˇrirozen´e ˇc´ıslo n. Program rozdˇel´ı interval ha, bina nstejn´ych d´ılk˚u. V dˇelic´ıch bodech xi,i= 0,1,...,n, spoˇcte f(xi) a vr´at´ı seznam subinterval˚u (xk, xk+1), v nichˇz plat´ı f(xk)f(xk+1)<0. (b) Naprogramujte metodu p˚ulen´ı intervalu. Vstupem bude funkce f(x), interval ha, bia poˇzadovan´a d´elka d=|bk−ak|v´ysledn´eho intervalu hak, bki. (c) Propojte (a) s (b); dostanete n´astroj pro hled´an´ı kvalitn´ıch poˇc´ateˇcn´ıch aproximac´ı. (d) Naprogramujte metodu seˇcen a teˇcen vˇcetnˇe testov´an´ı poˇzadovan´e pˇresnosti ε. Poˇc´ateˇcn´ı aproximaci urˇcete pomoc´ı (c) jako x0= (ak+bk)/2. 5.4. Jak´e v´ysledky pˇrinese program ze cviˇcen´ı 5.3 pro vstupn´ı data: (a) f(x) = x3+ 5x2−10, ha, bi=h−5,3i,n= 5, d= 0,1, ε= 10−6? (b) f(x) = x5+ 2x4−x3−2x2+ 0,1, ha, bi=h−3,3i,n= 5, d= 0,1, ε= 10−6? (c) f(x) = x5+ 2x4−x3−2x2+ 0,1, ha, bi=h−3,3i,n= 10, d= 0,1, ε= 10−6? [ (a) Najde vˇsechny koˇreny x∗ 1 . =−4,507903, x∗ 2 . =−1,755640, x∗ 3 . = 1,263543; (b) najde jen koˇreny x∗ 1 . =−2,008176, x∗ 2 . =−0,945472, x∗ 3 . = 0,982479; (c) najde vˇsechny koˇreny x∗ 1 . =−2,008176, x∗ 2 . =−0,945472, x∗ 3 . = 0,982479, x∗ 4 . =−0,246397, x∗ 5 . = 0,217566.] 5.5. Ovˇeˇrte, ˇze Newtonova metoda pro f(x) = (x−1)2konverguje k x∗= 1 jen line´arnˇe. [xk+1 −1 = 1 2(xk−1). ] 5.6. Hledejte koˇren rovnice f(x) := x2+ lnx−10/x = 0 na intervalu I=h1,4imetodou prost´e iterace. Zkoumejte na poˇc´ıtaˇci konvergenci pro iteraˇcn´ı funkce (a) g(x) = e10/x−x2, (b) g(x) = 10/(x2+ lnx), (c) g(x) = p10/x −lnx, (d) g(x) = 3 √10 −xlnx. Poˇc´ateˇcn´ı iteraci x0zvolte vˇzdy ve stˇredu intervalu, tj. x0= 2,5. Ovˇeˇrte, ˇze f(x) = 0 ⇐⇒ x=g(x) a odhadnˇete q= max x∈I|g′(x)|. Jsou v´ysledky v souladu s vˇetou o konvergenci? [ (a) nekonverguje, q. = 9,72 ·104; (b) nekonverguje, q. = 30; (c) konverguje, q. = 1,74; (d) konverguje rychleji, q. = 0,57. ] 5.7. Newtonovou metodou ˇreˇste soustavy rovnic na 6 desetinn´ych cifer pˇresnˇe. (a) 2 cos(xy) = 1 2 sin(x+y) = 1 (b) 2 cos(xy)−sin x= 0 2xsin y−3ysin x=−1(c) x2y−exy =−2 ln(x+ 1)y−y2/x = 0 Jako starovac´ı hodnoty pouˇzijte (a) x0= 2,5, y0= 0,25; (b) x0= 1, y0= 1; (c) x0= 1, y0= 1. Bl´ızkost koˇrene ovˇeˇrte graficky. [ (a) x∗. = 2,125254, y∗. = 0,492740; (b) x∗. = 1,023402, y∗. = 1,103856; (c) x∗. = 1,256558, y∗. = 1,022638.] 96 6. Optimalizace Optimalizaˇcn´ı ´ulohy se zab´yvaj´ı v´ybˇerem nejlepˇs´ıch ˇreˇsen´ı z dan´e mnoˇziny moˇzn´ych ˇreˇsen´ı. Matematicky m˚uˇzeme optimalizaˇcn´ı ´ulohu formulovat jako nalezen´ı prvku x∗∈M takov´eho, ˇze pro libovoln´y prvek x∈Mplat´ı f(x∗)≤f(x)∀x∈M , (6.1) kde f:M7→ Rje minimalizovan´a (nˇekdy se tak´e ˇr´ık´a ´uˇcelov´a nebo c´ılov´a nebo kriteri´aln´ı) funkce aMje mnoˇzina pˇr´ıpustn´ych ˇreˇsen´ı. Jestliˇze pˇr´ıpustn´ym ˇreˇsen´ım m˚uˇze b´yt kaˇzd´y bod x= (x1, x2,...,xn)Tn-rozmˇern´eho Euklidova prostoru Rn, tj. M=Rn, hovoˇr´ıme onepodm´ınˇen´e optimalizaci. O funkci fbudeme pˇredpokl´adat, ˇze je spojit´a (pˇr´ıpadnˇe i se sv´ymi prvn´ımi a dalˇs´ımi derivacemi). Optimalizaˇcn´ı ´uloha (6.1) se naz´yv´a ´ulohou glob´aln´ı optimalizace. My se v t´eto kapitole omez´ıme na jednoduˇsˇs´ı ´ulohu lok´aln´ı optimalizace spoˇc´ıvaj´ıc´ı v nalezen´ı lok´aln´ıho minima, tj. prvku x∗∈Mtakov´eho, ˇze plat´ı f(x∗)≤f(x)∀x∈M∩O(x∗),(6.2) kde O(x∗) je nˇejak´e okol´ı bodu x∗. Podrobnˇeji si vˇsimneme dvou speci´aln´ıch ´uloh: v odstavci 6.1 se sezn´am´ıme s metodami pro minimalizaci funkce jedn´e promˇenn´e na intervalu ha, bia v odstavci 6.2 se budeme vˇenovat metod´am nepodm´ınˇen´e minimalizace funkce v´ıce promˇenn´ych. Pozn´amka. Urˇcen´ı maxima funkce g(x) m˚uˇzeme pˇrev´est na ´ulohu urˇcen´ı minima funkce f(x) = −g(x).  6.1. Jednorozmˇern´a minimalizace V tomto odstavci uvedeme metody pro pˇribliˇzn´e urˇcen´ı bodu x∗lok´aln´ıho minima funkce f(x) na intervalu ha, bi. Jestliˇze je funkce fna intervalu ha, biunimod´aln´ı, tj. kdyˇz m´a v intervalu ha, bijedin´e minimum, pak ho (d´ale uveden´ymi metodami pˇribliˇznˇe) najdeme. Pokud vˇsak m´a funkce fna intervalu ha, biv´ıce lok´aln´ıch minim, najdeme jedno z nich. Intervalovou bisekci pouˇz´ıt nem˚uˇzeme: i kdyˇz zn´ame f(a), f(b) a f((a+b)/2), nedok´aˇzeme rozhodnout, ve kter´e polovinˇe intervalu ha, biminimum leˇz´ı. Pouˇz´ıt lze intervalovou trisekci. Necht ’h= (b−a)/3, takˇze u=a+hav=b−hdˇel´ı interval na tˇri stejn´e ˇc´asti. Pˇredpokl´adejme, ˇze f(u)< f(v). Pak minimum jistˇe leˇz´ı vlevo od v, takˇze bnahrad´ıme pomoc´ı v. T´ım se d´elka intervalu (obsahuj´ıc´ıho minimum) zkr´at´ı na dvˇe tˇretiny sv´e p˚uvodn´ı d´elky. Bod use vˇsak stane stˇredem nov´eho intervalu a nebude proto v dalˇs´ım kroku vyuˇziteln´y. Funkci ftedy mus´ıme vyhodnocovat v kaˇzd´em kroku dvakr´at. To je neefektivn´ı. Metoda zlat´eho ˇrezu je zaloˇzena na ˇsikovnˇejˇs´ım v´ybˇeru dˇelic´ıch bod˚u uav. Necht ’ h=(b−a), kde je ˇc´ıslo o nˇeco vˇetˇs´ı neˇz 1/3, jehoˇz pˇresnou hodnotu teprve urˇc´ıme. Pak body u=a+hav=b−hdˇel´ı interval ha, bina tˇri nestejn´e ˇc´asti. V prvn´ım kroku vyhodnot´ıme f(u) a f(v). Pˇredpokl´adejme, ˇze f(u)< f(v). Pak v´ıme, ˇze minimum je 97 mezi aav. Nahrad´ıme bpomoc´ı va proces opakujeme. Kdyˇz zvol´ıme spr´avnou hodnotu , bod ubude ve spr´avn´e pozici pouˇziteln´e v pˇr´ıˇst´ım kroku. Po prvn´ım kroku se tak funkce fbude vyhodnocovat v kaˇzd´em kroku uˇz jen jednou. Jak tedy zvolit ? Tak, aby bod uhr´al v redukovan´em intervalu ha, vistejnou roli jako bod vv p˚uvodn´ım intervalu ha, bi, tj. aby pomˇer d´elky intervalu ha, uik d´elce intervalu ha, vibyl stejn´y jako pomˇer d´elky intervalu ha, vik d´elce intervalu ha, bi, u−a v−a=v−a b−a⇐⇒  1−=1− 1⇐⇒ 2−3+ 1 = 0 . Vyhovuj´ıc´ı ˇreˇsen´ı je = (3 −√5)/2≈0,382,kde τ= 1 −= (√5−1)/2≈0,618 je ˇc´ıslo zn´am´e jako pomˇer zlat´eho ˇrezu. a bu ρ(b−a) v ρ(b−a) (1−ρ)(b−a) (1−ρ)(b−a) Obr. 6.1: Zlat´y ˇrez Pozn´amka (O zlat´em ˇrezu). ˇ R´ık´ame, ˇze bod dˇel´ı interval v pomˇeru zlat´eho ˇrezu, kdyˇz dva novˇe vznikl´e subintervaly maj´ı tuto vlastnost: pomˇer d´elky kratˇs´ıho subintervalu k d´elce delˇs´ıho subintervalu je stejn´y jako pomˇer d´elky delˇs´ıho subintervalu k d´elce cel´eho intervalu. Z v´yˇse uveden´e konstrukce je zˇrejm´e, ˇze bod u(ale tak´e bod v) dˇel´ı interval ha, biv pomˇeru zlat´eho ˇrezu. Pˇripomeˇnme si, ˇze s ˇc´ıslem τjsme se setkali jiˇz v kapitole 5.2, kde jsme uvedli, ˇze rychlost konvergence metody seˇcen p= 1 + τ. = 1,618.  Zat´ım jsme pˇredpokl´adali, ˇze f(u)< f(v). V opaˇcn´em pˇr´ıpadˇe, tj. kdyˇz f(u)≥f(v), leˇz´ı minimum v intervalu hu, bi, takˇze anahrad´ıme pomoc´ı u. Snadno ovˇeˇr´ıme, ˇze v redukovan´em intervalu hu, bibude m´ıt vstejnou roli jako mˇelo uv p˚uvodn´ım intervalu ha, bi, takˇze hodnotu funkce fna redukovan´em intervalu opˇet staˇc´ı poˇc´ıtat jen jednou. D´elka redukovan´eho intervalu je τ-kr´at menˇs´ı neˇz d´elka p˚uvodn´ıho intervalu. Z v´ychoz´ıho intervaluha0, b0i=ha, bitak postupnˇe sestroj´ıme intervaly ha1, b1i ⊃ ha2, b2i ⊃ ..., kter´e obsahuj´ı minimum a jejichˇz d´elka je v kaˇzd´em kroku redukov´ana faktorem τ. Na v´ychoz´ım intervalu ha0, b0iurˇc´ıme u0=a+(b−a), v0=b−(b−a) a vypoˇcteme f(u0), f(v0). Interval hak+1, bk+1i,k= 0,1,..., dostaneme pomoc´ı ak,bk,uk,vka jiˇz dˇr´ıve vypoˇcten´ych hodnot f(uk), f(vk) takto: 1) kdyˇz f(uk)< f(vk), pak ak+1 := ak,bk+1 := vk,vk+1 := uk,uk+1 := ak+1 +bk+1 −vk+1 a vypoˇcteme f(uk+1); 2) v opaˇcn´em pˇr´ıpadˇe, tj. kdyˇz f(uk)≥f(vk), provedeme ak+1 := uk,bk+1 := bk,uk+1 := vk,vk+1 := ak+1 +bk+1 −uk+1 98 a vypoˇcteme f(vk+1). akbk ukvk ak+1 bk+1 akbk ukvk ak+1 bk+1 Obr. 6.2: Metoda zlat´eho ˇrezu Po kkroc´ıch leˇz´ı minimum v intervalu Ik:= hak, bkid´elky |Ik|=bk−ak=τ(bk−1−ak−1) = ···=τk(b0−a0). Stˇred xk+1 intervalu hak, bkiaproximuje minimum x∗s chybou |xk+1 −x∗| ≤ 1 2(bk−ak) = 1 2τk(b0−a0).(6.3) Pro k→ ∞ zˇrejmˇe |Ik| → 0 a xk→x∗. Konvergence metody zlat´eho ˇrezu je pomˇernˇe pomal´a. Proto je v bl´ızkosti minima ´uˇceln´e pˇrej´ıt na rychleji konvergentn´ı metodu. Jednou z moˇznost´ı je metoda kvadratick´e interpolace, s n´ıˇz se tak´e sezn´am´ıme. Pˇredt´ım ale Pˇr´ıklad 6.1. Metodou zlat´eho ˇrezu urˇc´ıme minimum funkce f(x) = x4−3x3+x+7. Jako poˇc´ateˇcn´ı zvol´ıme interval ha0, b0i=h1,3i. V´ypoˇcet provedeme s pˇresnost´ı ε= 10−3, tj. kdyˇz bk−ak<2ε, poloˇz´ıme xk+1 = (ak+bk)/2. V´ypoˇcet je zaznamen´an v n´asleduj´ıc´ı tabulce. Podtrˇzen´ım jsou vyznaˇceny ty vnitˇrn´ı body, v nichˇz se poˇc´ıt´a hodnota ´uˇcelov´e funkce. k akukvkbkf(uk)f(vk) 0 1,0000 1,7639 2,2361 3,0000 1,979900 ≥0,695048 1 1,7639 2,2361 2,5279 3,0000 0,695048 <1,901312 2 1,7639 2,0557 2,2361 2,5279 0,852324 ≥0,695048 3 2,0557 2,2361 2,3475 2,5279 0,695048 <0,906510 4 2,0557 2,1672 2,2361 2,3475 0,690296 <0,695048 5 2,0557 2,1246 2,1672 2,2361 0,729249 ≥0,690296 . . . 14 2,1973 2,1982 2,1988 2,1997 0,681572 <0,681575 15 2,1973 2,1988 Poˇzadovan´a pˇresnost je dosaˇzena pro k= 15, takˇze (po zaokrouhlen´ı na 3 desetinn´e cifry) x16 . = 2,198. Hodnota ´uˇcelov´e funkce se poˇc´ıt´a celkem 16-kr´at. Protoˇze pˇresn´a hodnota x∗. = 2,198266, m´a x16 := 2,198 vˇsechny cifry platn´e.  99 Metoda kvadratick´e interpolace. Pˇredpokl´adejme, ˇze minimum leˇz´ı v intervalu hak, bki, a ˇze v nˇejak´em jeho vnitˇrn´ım bodˇe ckhodnota funkce fnepˇres´ahne hodnoty f(ak), f(bk) v krajn´ıch bodech ak,bk, tj. ˇze pro ak< ck< bkplat´ı f(ck)≤min{f(ak); f(bk)}.(6.4) Body [ak, f(ak)], [ck, f(ck)] a [bk, f(bk)] proloˇz´ıme parabolu P2(x) (kvadratick´y interpolaˇcn´ı polynom) a bod xk+1 jej´ıho minima povaˇzujeme za dalˇs´ı aproximaci x∗. Vzorec pro v´ypoˇcet xk+1 dostaneme ˇreˇsen´ım line´arn´ı rovnice P′ 2(xk+1) = 0. D´a se uk´azat, ˇze xk+1 =ck−1 2 (ck−ak)2[f(ck)−f(bk)] −(ck−bk)2[f(ck)−f(ak)] (ck−ak)[f(ck)−f(bk)] −(ck−bk)[f(ck)−f(ak)] .(6.5) Z (6.4) plyne, ˇze xk+1 ∈(ak, bk). Kdyˇz n´ahodou xk+1 =ck, vloˇz´ıme do xk+1 jin´y vnitˇrn´ı bod intervalu hak, bki. S t´ımto slab´ym m´ıstem metody kvadratick´e interpolace (a s dalˇs´ımi, zde nezm´ınˇen´ymi nedostatky) se ´uspˇeˇsnˇe vyrovn´av´a Brentova metoda. Struˇcn´a zm´ınka o n´ı je uvedena v n´asleduj´ıc´ım textu. Z bod˚u ak,bk,ckaxk+1 pak vybereme nov´y interval (ak+1, bk+1) obsahuj´ıc´ı minimum a bod ck+1 splˇnuj´ıc´ı podm´ınku (6.4), tentokr´at pro index k+ 1. Postupujeme podle n´asleduj´ıc´ıch pravidel: A : xk+1 < ckaf(xk+1)< f(ck) =⇒ak+1 =ak,ck+1 =xk+1,bk+1 =ck, B : xk+1 < ckaf(xk+1)≥f(ck) =⇒ak+1 =xk+1,ck+1 =ck,bk+1 =bk, C : ck< xk+1 af(ck)< f(xk+1) =⇒ak+1 =ak,ck+1 =ck,bk+1 =xk+1, D : ck< xk+1 af(ck)≥f(xk+1) =⇒ak+1 =ck,ck+1 =xk+1,bk+1 =bk. akbk xk+1 ck ak+1 ck+1 bk+1 A akbk xk+1 ck ak+1 ck+1 bk+1 B akbk ckxk+1 ak+1 ck+1 bk+1 C akbk ckxk+1 ak+1 ck+1 bk+1 D Obr. 6.3: Metoda kvadratick´e interpolace V´ypoˇcet ukonˇc´ıme a xk+1 povaˇzujeme za dostateˇcnˇe dobrou aproximaci minima x∗, kdyˇz je splnˇeno nˇekter´e ze stop kriteri´ı |xk+1 −xk|< ε , |xk+1 −xk|< ε|xk|,|f(xk)−f(xk+1)|< ε , |f(xk)−f(xk+1)|< ε|f(xk)|, 100 kde εje pˇredepsan´a tolerance. V´ypoˇcet xk+1,k= 0,1,..., vyˇzaduje k+ 3 vyhodnocen´ı ´uˇcelov´e funkce: f(a0), f(c0), f(b0) a d´ale f(xi), i= 1,2,...,k. Pokud metoda kvadratick´e interpolace konverguje, pak je rychlost jej´ı konvergence superline´arn´ı ˇr´adu p≈1,324, viz [7]. Pˇr´ıklad 6.2. Minimum funkce f(x) = x4−3x3+x+7 urˇc´ıme metodou kvadratick´e interpolace. V´ypoˇcet ukonˇc´ıme, kdyˇz |xk+1 −xk|<10−5. V´ypoˇcet je zaznamen´an v n´asleduj´ıc´ı tabulce. V posledn´ım sloupci je uvedeno, kter´y z pˇr´ıpad˚u A,B,C,D nast´av´a, a podtrˇzen´ım a0= 2,000000 x1= 2,204918 < c0= 2,500000 b0= 3,000000 f(x1)< f(c0) =⇒A a1= 2,000000 x2= 2,180689 < c1= 2,204918 b1= 2,500000 f(x2)≥f(c1) =⇒B a2= 2,180689 x3= 2,197322 < c2= 2,204918 b2= 2,500000 f(x3)< f(c2) =⇒A a3= 2,180689 c3= 2,197322 < x4= 2,198232 b3= 2,204918 f(c3)≥f(x4) =⇒D a4= 2,197322 c4= 2,198232 < x5= 2,198264 b4= 2,204918 f(c4)≥f(x5) =⇒D a5= 2,198232 c5= 2,198264 < x6= 2,198265 b5= 2,204918 jsou (postupnˇe zleva doprava) vyznaˇceny body ak+1 < ck+1 < bk+1. Z tabulky je zˇrejm´e, ˇze poˇzadovan´a pˇresnost byla dosaˇzena pro x6. = 2,19827 (vˇsechny cifry jsou platn´e). Hodnota ´uˇcelov´e funkce se poˇc´ıtala celkem 8-kr´at.  Brentova metoda je kombinovan´a metoda, kter´a v sobˇe spojuje spolehlivost metody zlat´eho ˇrezu a rychlou konvergenci metody kvadratick´e interpolace. Popis Brentovy metody lze naj´ıt napˇr. v [18], viz funkce brent, nebo v [15], viz funkce fmintx. Brent˚uv algoritmus je tak´e z´akladem funkce fminbnd pro jednorozmˇernou minimalizaci v MATLABu. 6.2. Minimalizace funkce v´ıce promˇenn´ych Nelderova-Meadova metoda zn´am´a tak´e jako metoda simplex˚u je popul´arn´ı metoda nepouˇz´ıvaj´ıc´ı derivace ´uˇcelov´e funkce. Patˇr´ı mezi tzv. komparativn´ı metody, coˇz jsou metody, kter´e hledaj´ı minimum ´uˇcelov´e funkce fporovn´av´an´ım jej´ıch hodnot v urˇcit´ych vybran´ych bodech prostoru Rn. V pˇr´ıpadˇe metody simplex˚u jsou vybran´ymi body vrcholy simplexu (pro n= 2 troj´uheln´ıka, pro n= 3 ˇctyˇrstˇenu). Hlavn´ı myˇslenka jednoho kroku metody je jednoduch´a: mezi vrcholy x0,x1,...,xnsimplexu vybereme nejhorˇs´ı vrchol xw(v angliˇctinˇe worst), v nˇemˇz ´uˇcelov´a funkce nab´yv´a nejvˇetˇs´ı hodnotu, a nahrad´ıme ho lepˇs´ım vrcholem ˆ x, v nˇemˇz je hodnota ´uˇcelov´e funkce menˇs´ı. Vrchol ˆ xhled´ame na polopˇr´ımce, kter´a vych´az´ı z nejhorˇs´ıho vrcholu xwa proch´az´ı tˇeˇziˇstˇem ¯ xzb´yvaj´ıc´ıch vrchol˚u. Nejlepˇs´ı z nich oznaˇc´ıme xb(v angliˇctinˇe best), tj. xbje ten z vrchol˚u x0,x1,...,xn, v nˇemˇz ´uˇcelov´a funkce nab´yv´a nejmenˇs´ı hodnotu. Prvn´ı pokus, jak vybrat ˆ x, oznaˇcujeme jako reflexi: bod xr=¯ x+ (¯ x−xw), je obraz bodu xwve stˇredov´e soumˇernosti se stˇredem ¯ x. Kdyˇz je f(xr)< f(xb), pak to znamen´a, ˇze pokles hodnot na polopˇr´ımce xw¯ xje znaˇcn´y, a proto zkus´ıme postoupit po t´eto polopˇr´ımce jeˇstˇe d´al, do bodu xe=¯ x+2(¯ x−xw). V´ybˇer bodu xeb´yv´a oznaˇcov´an jako expanze. Kdyˇz f(xe)< f(xb), pak ˆ x=xe, tj. nejhorˇs´ı vrchol xwnahrad´ıme bodem xe. 101 Jestliˇze f(xe)≥f(xb), zkus´ıme pouˇz´ıt alespoˇn bod xr. Podm´ınkou pro jeho zaˇrazen´ı do simplexu je splnˇen´ı podm´ınky f(xr)< f(xg) pro nˇekter´y vrchol xgjin´y neˇz nejhorˇs´ı, tj. pro xg6=xw(index gpˇripom´ın´a anglick´e sl˚uvko good). Pokud takov´a podm´ınka plat´ı, bereme ˆ x=xrm´ısto p˚uvodn´ıho xw. 1 2 3 4 5 6 7 8 xb xg xw xr xe xce xci xb xg xw 0) 1) 2) 3) 4) 5) Obr. 6.4: Nelderova-Meadova metoda: 0) origin´aln´ı troj´uheln´ık, 1) expanze, 2) reflexe, 3) vnˇejˇs´ı kontrakce, 4) vnitˇrn´ı kontrakce, 5) redukce Kdyˇz nevyhovuje xeani xr, zkus´ıme naj´ıt bod ˆ xna ´useˇcce s koncov´ymi body xw,xr tak, aby f(ˆ x)<min{f(xw); f(xr)}. Konkr´etnˇe postupujeme takto: a) pokud f(xr)< f(xw), zkus´ıme bod xce =1 2(¯ x+xr) (leˇz´ı bl´ıˇze k bodu xr), a kdyˇz f(xce)< f(xr), pak ˆ x=xce, takˇze provedeme xw:= xce. b) jestliˇze f(xr)≥f(xw), zkus´ıme bod xci =1 2(¯ x+xw) (leˇz´ı bl´ıˇze k bodu xw), a pokud f(xci)< f(xw), pak ˆ x=xci, tj. provedeme xw:= xci. V´ybˇer bodu xce resp. xci oznaˇcujeme jako kontrakci (v doln´ım indexu: p´ısmeno c pˇripom´ın´a anglick´e slovo contraction, p´ısmeno eanglick´e slovo external (xce leˇz´ı vnˇe p˚uvodn´ıho simplexu) a p´ısmeno ipˇripom´ın´a anglick´e slovo internal (xci leˇz´ı uvnitˇr p˚uvodn´ıho simplexu)). Kdyˇz nevyhovuje xe,xr,xce ani xci, usoud´ıme, ˇze vrchol xb, v nˇemˇz ´uˇcelov´a funkce nab´yv´a sv´e nejmenˇs´ı hodnoty, je bl´ızko minima. Proto provedeme redukci simplexu: vrchol 102 xbv simplexu z˚ustane a zb´yvaj´ıc´ı vrcholy xise posunou do stˇred˚u ´useˇcek xbxi. Simplex se tedy st´ahne k nejlepˇs´ımu vrcholu xb. Transformaci simplexu, pˇredstavuj´ıc´ı jeden krok Nelderovy-Meadovy metody, pop´ıˇseme v pˇeti bodech takto: 1) expanze : f(xr)< f(xb) a nav´ıc f(xe)< f(xb) =⇒xw:= xe 2) reflexe : f(xr)< f(xg) pro nˇejak´y vrchol xg6=xw=⇒xw:= xr 3) vnˇejˇs´ı kontrakce : f(xr)< f(xw) a nav´ıc f(xce)< f(xr) =⇒xw:= xce 4) vnitˇrn´ı kontrakce : f(xr)≥f(xw) a nav´ıc f(xci)< f(xw) =⇒xw:= xci 5) redukce : xi:= 1 2(xb+xi) pro vˇsechna xi6=xb Body 1 aˇz 5 proch´az´ıme postupnˇe shora dol˚u. Kdyˇz nˇekter´a z podm´ınek v bodech 1 aˇz 4 nen´ı splnˇena, pˇrejdeme na n´asleduj´ıc´ı bod. Kdyˇz splnˇena je, provedeme n´ahradu xwpodle pˇr´ıkazu za ˇsipkou a transformace je hotova. Nen´ı-li splnˇena podm´ınka v ˇz´adn´em z bod˚u 1 aˇz 4, provedeme redukci podle bodu 5. Na zaˇc´atku v´ypoˇctu je d´ana poˇc´ateˇcn´ı aproximace x0= (x(0) 1, x(0) 2,...,x(0) n)Ta mal´e ˇc´ıslo δ. Dalˇs´ı vrcholy xistartovac´ıho simplexu odvod´ıme z vrcholu x0tak, ˇze k jeho i-t´e sloˇzce x(0) ipˇriˇcteme ˇc´ıslo δ, tj. xi= (x(0) 1,...,x(0) i+δ,...,x(0) n)T,i= 1,2,...,n. Simplex opakovanˇe transformujeme. V´ypoˇcet ukonˇc´ıme a vrchol xbpovaˇzujeme za dostateˇcnˇe dobrou aproximaci minima x∗, pokud jsou vrcholy simplexu navz´ajem dosti bl´ızko a funkˇcn´ı hodnoty v nich se m´alo liˇs´ı, tj. kdyˇz pro zadan´e tolerance ε1,ε2plat´ı kxi−xbk< ε1a souˇcasnˇe |f(xi)−f(xb)|< ε2pro kaˇzd´y vrchol xi6=xb.(6.6) Nelderova-Meadova metoda je heuristick´a metoda (heuristick´y nebo-li objevovac´ı je takov´y postup, kter´y je zaloˇzen nejenom na logick´em uvaˇzov´an´ı a zkuˇsenostech, ale tak´e na pozorov´an´ı a experimentov´an´ı). Metoda je vhodn´a pro minimalizaci funkc´ı menˇs´ıho poˇctu promˇenn´ych, ˇreknˇeme pro n≤10. Pˇrestoˇze o konvergenci metody je toho zn´amo jen velmi m´alo, praxe hovoˇr´ı v jej´ı prospˇech: metoda je aˇz pˇrekvapivˇe ´uspˇeˇsn´a. Proto je povaˇzov´ana za metodu spolehlivou nebo-li robustn´ı. Podstatnou nev´yhodou metody je to, ˇze je pomal´a, zejm´ena v bl´ızkosti minima. Dalˇs´ı minus pˇredstavuje velk´y objem v´ypoˇct˚u. Pˇresto nejde o mrtvou metodu, coˇz je zˇrejm´e napˇr. z toho, ˇze je implementov´ana v MATLABu jako funkce fminsearch. Pˇr´ıklad 6.3. Funkce f(x, y) = 70 [(x−2)4+(x−2y)2] nab´yv´a minima pro x∗= 2, y∗= 1. V´ypoˇcet provedeme metodou Neldera-Meada, v n´ıˇz zvol´ıme x0= (2,1; 0,7)Taδ= 0,1. V´ypoˇcty jsou prov´adˇeny pˇresnˇe, do tabulky 6.1 vˇsak (kv˚uli ´uspoˇre m´ısta) zapisujeme hodnoty zaokrouhlen´e na dvˇe desetinn´a m´ısta. Pro kaˇzd´y bod zapisujeme ve sloupci pod sebou x-ovou souˇradnici, y-ovou souˇradnici a funkˇcn´ı hodnotu. V r´ameˇcku uv´ad´ıme novˇe urˇcen´y lepˇs´ı bod ˆ x(a funkˇcn´ı hodnotu v nˇem). Pro ε1= 10−4aε2= 10−8v´ypoˇcet konˇc´ı po 47 kroc´ıch a xb. = (1,999953; 0,999976)T. Metoda nejvˇetˇs´ıho sp´adu je z´akladn´ı minimalizaˇcn´ı metoda, kter´a pouˇz´ıv´a derivace ´uˇcelov´e funkce. Takov´e metody se naz´yvaj´ı gradientn´ı. Zaˇcneme t´ım, ˇze si vysvˇetl´ıme obecn´y princip sp´adov´e metody. Pˇredpokl´adejme tedy, ˇze jsme v bodu xka chceme se dostat bl´ıˇze k minimu. Zvol´ıme smˇer dk, v nˇemˇz funkce f 103 krok typ xwxgxb¯ x xrxexce xci 1 expanze 2,20 0,70 44,91 2,10 0,70 34,31 2,10 0,80 17,51 2,10 0,75 − 2,00 0,80 11,20 1,90 0,85 2,81 2 reflexe 2,10 0,70 34,31 2,10 0,80 17,51 1,90 0,85 2,81 2,00 0,82 − 1,90 0,95 0,01 1,80 1,08 8,69 3vnˇejˇs´ı kontrakce 2,10 0,80 17,51 1,90 0,85 2,81 1,90 0,95 0,01 1,90 0,90 − 1,70 1,00 6,87 1,80 0,95 0,81 4vnitˇrn´ı kontrakce 1,90 0,85 2,81 1,80 0,95 0,81 1,90 0,95 0,01 1,85 0,95 − 1,80 1,05 6,41 1,88 0,90 0,41 Tab. 6.1: Pˇr´ıklad 6.3 kles´a, a na polopˇr´ımce xk+λdk,λ≥0, vybereme bod xk+1 =xk+λkdk,(6.7) v nˇemˇz f(xk+1)< f(xk). Smˇerov´y vektor dknaz´yv´ame sp´adov´y (ve smˇeru dkhodnota ´uˇcelov´e funkce fpad´a dol˚u), odtud sp´adov´a metoda. ˇ C´ıslo λkse naz´yv´a parametr d´elky kroku (je-li dkjednotkov´y vektor, tj. kdyˇz kdkk2= 1, pak λk=kxk+1 −xkk2je vzd´alenost bod˚u xkaxk+1, tj. λkje d´elka kroku). λkdostaneme minimalizac´ı funkce ϕ(λ) = f(xk+ λdk) pro λ≥0. Minimum λkfunkce ϕ(λ) urˇc´ıme pˇribliˇznˇe pomoc´ı nˇekolika m´alo krok˚u vhodn´e metody jednorozmˇern´e minimalizace (pˇresn´a minimalizace je zbyteˇcn´y pˇrepych, nebot ’prostˇrednictv´ım λkurˇcujeme jen meziv´ysledek xk+1 na cestˇe k minimu x∗). Urˇcen´ı λktedy vyj´adˇr´ıme z´apisem λk. = argmin λ>0 ϕ(λ),kde ϕ(λ) = f(xk+λdk).(6.8) Nyn´ı se vˇenujme uˇz vlastn´ı metodˇe nejvˇetˇs´ıho sp´adu. Je zn´amo, ˇze funkce f(x) nejrychleji kles´a ve smˇeru z´aporn´eho gradientu. Oznaˇc´ıme-li tedy gradient jako g(x)≡ ∇f(x) = ∂f(x) ∂x1 ,∂f(x) ∂x2 ,...,∂f(x) ∂xnT , pak v metodˇe nejvˇetˇs´ıho sp´adu vol´ıme jako sp´adov´y vektor dk=−g(xk).(6.9) V´ypoˇcet ukonˇc´ıme a xk+1 povaˇzujeme za uspokojivou aproximaci minima x∗, kdyˇz kg(xk+1)k< ε , (6.10) 104 x3 x2 x1 x0 Obr. 6.5: Princip sp´adov´ych metod kde εje pˇredepsan´a tolerance (pˇripomeˇnme, ˇze v minimu g(x∗) = o). V poˇc´ateˇcn´ı f´azi v´ypoˇctu, kdyˇz jsme od minima jeˇstˇe dosti daleko, doch´az´ı obvykle k pomˇernˇe rychl´emu poklesu hodnot ´uˇcelov´e funkce. Zato v bl´ızkosti minima je konvergence pomal´a, jen line´arn´ı, tj. plat´ı kxk+1 −x∗k ≤ Ckxk−x∗k, kde konstanta Cje sice menˇs´ı neˇz jedna, ale ˇcasto jen nepatrnˇe. Cik-cak efekt. Pokud λkvypoˇcteme pˇresnˇe, m´a funkce ϕ(λ) v bodˇe λkminimum, a proto 0 = ϕ′(λk) = n X i=1 ∂f(xk+1) ∂xi d(k) i=− n X i=1 d(k+1) id(k) i=−dT k+1dk, tj. smˇerov´e vektory dk+1 adkjsou navz´ajem kolm´e. Uk´aˇzeme si, ˇze pr´avˇe tato vlastnost m˚uˇze b´yt pˇr´ıˇcinou velmi pomal´e konvergence. Pro jednoduchost pˇredpokl´adejme, ˇze minimalizujeme funkci dvou promˇenn´ych. V tom pˇr´ıpadˇe si m˚uˇzeme vypomoci jednoduchou pˇredstavou: nach´az´ıme se v ter´enu a chceme naj´ıt nejniˇzˇs´ı bod, tj. dno nˇejak´e prohlubnˇe. Pˇredpokl´adejme, ˇze prohlubeˇn m´a tvar podlouhl´e zahnut´e rokle. Metoda nejvˇetˇs´ıho sp´adu n´as nasmˇeruje z poˇc´ateˇcn´ıho stanoviˇstˇe x0dol˚u kolmo k vrstevnici f(x) = f(x0). Sestupujeme tak dlouho, dokud ter´en kles´a. V nejniˇzˇs´ım m´ıstˇe je dalˇs´ı stanoviˇstˇe x1. Zde se zastav´ıme, otoˇc´ıme se o 90 stupˇn˚u a sestupujeme znovu dol˚u (kolmo k vrstevnici f(x) = f(x1)) do dalˇs´ıho stanoviˇstˇe x2 atd. Je zˇrejm´e, ˇze vzd´alenosti mezi jednotliv´ymi stanoviˇsti se budou postupnˇe zkracovat. V bl´ızkosti minima bude naˇse putov´an´ı p˚usobit smˇeˇsnˇe: m´ısto abychom do nˇej doˇsli nˇekolika m´alo kroky, budeme se k nˇemu sp´ıˇse pl´ıˇzit neˇz bl´ıˇzit po trase tvoˇren´e ˇc´ım d´al 105