PA081: Programování numerických ^ - o vypočtu 2. Nelineární rovnice o jedné neznámé Aleš Křenek jaro 2019 PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 1/52 Obecná formulace problému ► hledáme řešení rovnice F (x) = G(x) kde alespoň jedna z funkcí F, G není lineární ► víceméně všechny numerické metody jsou iterační ► začínáme s odhadem řešení xo ► v každém kroku odhad postupně zpřesňujeme ► výpočet končí dosažením kritéria zastavení ► dosáhli jsme dostatečně přesné aproximace řešení PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol Obecná formulace problému ► hledáme řešení rovnice F (x) = G(x) kde alespoň jedna z funkcí F, G není lineární ► víceméně všechny numerické metody jsou iterační ► začínáme s odhadem řešení xo ► v každém kroku odhad postupně zpřesňujeme ► výpočet končí dosažením kritéria zastavení ► dosáhli jsme dostatečně přesné aproximace řešení ► anebo ► rovnice nemá řešení ► použitá metoda pro danou rovnici nefunguje ► špatně jsme odhadli Xo ► dosáhli jsme falešného řešení ► žádná metoda není dokonale univerzální ► bez jisté analýzy vlastností rovnice se neobejdeme Železniční příklad Kolejnice délky 2d je ohnutá do oblouku tak, že její konce jsou ve vzálenosti 2a. Jaká je vzdálenost středu kolejnice od spojnice krajních bodů? ► označíme R poloměr oblouku, a úhel poloviny úseče, r hledanou vzdálenost ► platí rovnice d = Ra Rsina = a r = R(l-cosa) ► dosazením a jednoduchou úpravou d . — sin a = a a ► hledáme pevný bod funkce ► kandidát na metodu prosté iterace PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol Metoda prosté iterace ► rovnice ve tvaru x = f(x) ► řešením je pevný bod funkce / ► počítáme opakovaně Xí+i = f {x i) ► až k dosažení kritéria konvergence PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol Metoda prosté iterace ► ► ► rovnice ve tvaru x = f(x) řešením je pevný bod funkce / počítáme opakovaně Xi+i = f(xi) ► až k dosažení kritéria konvergence ► kdy a proč to funguje? - věta o pevném bodě Je-li pro K c [R funkce /: K — K kontrakce, tj. existuje q g (0,1) tak, že M x,y g K: \f(x) - f(y) \ < q\x - y\, potom existuje x* g K takové, že pro libovolné Xo g K je x* limitou posloupnosti Xi+i =f(Xi). ► idea důkazu Xi+i = \f(Xi) <\Xi-X PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol Metoda prosté iterace Pro železniční příklad PA081: Programování numerických výpočtů A. Křenek ► je zobrazení a »— % sin ex kontrakce? a ► postačující podmínka a V x e K: f'(x) < 1 tedy cos a < — a d ► | sin ex = a bude mít řešení v (arccos f, -f) Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 5/52 Metoda prosté iterace Pro železniční příklad ► je zobrazení a »— ^ sin ex kontrakce? ► postačující podmínka Mx e K\ fr(x) < 1 tedy cos a < % a ► ^ siná = a bude mít řešení v (arccos %\) ► implementačně velmi jednoduchá metoda ► zdánlivě velmi speciální případ ► řešený problém lze často transformovat, aby splňoval podmínku kontrakce ► rychlost konvergence záleží na vlastnostech f(x) Pohyb planet ► poloha planety na eliptické dráze o poloosách a, b x = acosE - e y = avl - e2 sin£ kde e = Jl- b2 a* E = cot + e sin£ ► nelze snadno řešit analyticky ► jeden z nej známějších motivačních příkladů pro numerické řešení rovnic ► viz např. „Kepleťs equation" na Wikipedii PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 6/52 Hledání kořenů ► pro rovnici F(x) = G(x) uvažujeme f(x) = F(x) - G(x) ► provedeme separaci kořenů f(x) ► definiční obor f(x) rozdělíme na intervaly [a*, bt\ ► v každém [a*, bt\ má funkce právě jeden kořen ► f(at)f(bi) < 0 ► / je na [aubi] spojitá ► nepovede-li se separace dokonale, musíme být připraveni na následky ► vybereme vhodnou metodu hledání kořene na [au b i] ► stanovíme konkrétní podmínku ukončení PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 7/52 Separace kořenů Spojitá funkce ► máme štěstí a / je spojitá, platí varianta věty o střední hodnotě Je-li / na [a, b] spojitá a f(a)f(b) < 0, existuje x e [a, b] tak, že f(x) = 0 a / b 1 m PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol Separace kořenů Nespojitá funkce ► není-li / spojitá, ale je ohraničená ► „kříží" osu x v bodě nespojitosti ► numericky nerozlišitelné od kořene PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol Separace kořenů Nespojitá funkce PA081: Programování numerických výpočtů A. Křenek ► není-li / spojitá, ale je ohraničená ► „kříží" osu x v bodě nespojitosti ► numericky nerozlišitelné od kořene ► není-li ani ohraničená ► např. x-c ► některé metody najdou „kořen" v c ► snadno identifikovatelný problém Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 9/52 Separace kořenů Patologické případy u u 1 1 1 e f c ^r 1 1 ^^^^ 1 X\ d a x-i x$ ŕ PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol Separace kořenů Praktický postup ► základní analýza vlastností / je nezastupitelná ► netušíme-li vůbec, zdali má kořeny, kde jsou, kolik jich je atd., je problém zpravidla někde jinde PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 11/52 Separace kořenů Praktický postup ► základní analýza vlastností / je nezastupitelná ► netušíme-li vůbec, zdali má kořeny, kde jsou, kolik jich je atd., je problém zpravidla někde jinde ► algoritmus „look inward" ► alespoň nahrubo známe interval, kde by kořen měl být ► rozdělíme na n menších a postupně prohledáváme ► v případě neúspěchu můžeme opakovat s větším n Separace kořenů Praktický postup ► základní analýza vlastností / je nezastupitelná ► netušíme-li vůbec, zdali má kořeny, kde jsou, kolik jich je atd., je problém zpravidla někde jinde ► algoritmus „look inward" ► alespoň nahrubo známe interval, kde by kořen měl být ► rozdělíme na n menších a postupně prohledáváme ► v případě neúspěchu můžeme opakovat s větším n ► algoritmus „look outward" ► poslední zoufalý pokus ► počáteční interval expandujeme exponenciálně na tu stranu, kde se funkce blíží k ose x ► zpravidla funguje pro funkce, které mají pro x — ±00 opačné zneménko Metoda půlení intervalů Princip ► začínáme se separovaným kořenem ► pro daný interval [a, b] platí f(a)f(b) < 0 b+a ► není-li přesně kořen, platí právě jedna nerovnost b + a f (a) < 0 b + a f(b)<0 ► interval rozpůlíme a pokračujeme rekurzivně ► metoda vždy dokonverguje ke kořeni nebo k singularitě ► je-li jich více, odhalí jen jeden PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 12/52 Metoda půlení intervalů Konvergence ► je-li požadovaná přesnost e, je třeba n = log2 b - a iteračních kroku ► lineární rychlost konvergence ► v každém kroku přibývá konstatně platných číslic přesnosti ► kritéria ukončení ► jsou-li a, b řádově srovnatelná, nemá smysl více iterací než než počet bitů mantisy ► v opačném případě opět musíme vědět, co chceme ► standardně se používá zastavení při velikosti intervalu \a\ + \b\ kde e je přesnost daného datového typu ► metoda je robustní ale relativně pomalá PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 13/52 Metoda půlení intervalů PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 14/52 Newtonova metoda ► také Newton-Raphsonova metoda ► odvozena z Taylorova rozvoje x* = Xi + ô i f(xt + ô t) = f(Xi) + f'ixdôi + f iXi) 1 + 2! ► zanedbáme vyšší derivace j * ~i t jŕ X PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol Newtonova metoda Konvergence ► označíme x* kořen, a e j odchylky Xi - x* X* + Ci+i = Xi+i = Xi - *}Xl\ = X* + €i - 1 f'(Xi) f(Xj) 1 f'(xú tedy ei+i =ei-^-) ► Taylorův rozvoj e2J"(li) 0 = f(x*)=f(xi)-eif'(xi) + 2! kde §í g [0, C{] ► podělením f'(xi) a dosazením dostaneme ei+l = _eí 2/'(Xi) ► kvadratická konvergence ► v každém kroku se zdvojnásobí počet platných číslic PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 16/52 Newtonova metoda Nešvary ► špatná globální konvergence ► velká citlivost na lokální vlastnosti / mimo kořen 2 / > * t ► prakticky nepoužitelná sama o sobě ► nutno kombinovat s jinými metodami ► vhodná k „vyleštění" nahrubo nalezeného kořene 4 1 # M * w + m t X* *2 M m * M -ŕ m * # * f * M -ŕ w 4> PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 17/52 Newtonova metoda PA081: Programování numerických výpočtů A. Křenek float rtnewt(void (*funcd) (float, float *, float *) , float xl, float x2, float xacc) { int j; float df,dx,f,rtn; rtn=0.5*(xl+x2); for (j=l;j<=JMAX;j++) { (*funcd)(rtn,&f,&df); dx=f/df; rtn -= dx; if ((xl-rtn)*(rtn-x2) < 0.0) { /* chyba, utekli jsme */ return 0; } if (fabs(dx) < xacc) return rtn; } /* Přiliš mnoho iteraci */ return 0; Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 18/52 Halleyova metoda ► máme k dispozici i 2. derivaci ► zanedbáme až kubický člen Taylorova rozvoje f(Xí) f'(Xi)[l- .f(Xi).f"(Xi) 2/'(xí)2 PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 19/52 Halleyova metoda ► máme k dispozici i 2. derivaci ► zanedbáme až kubický člen Taylorova rozvoje __f(Xj) ► rychlost konvergence 3 ► za cenu výpočtu f"\ ► proč tedy nespočítat místo toho další krok N. metody (rychlost konvergence 4)? ► v mnoha případech lze recyklovat části výpočtu, /" je téměř zadarmo Halleyova metoda ► máme k dispozici i 2. derivaci ► zanedbáme až kubický člen Taylorova rozvoje f(Xi) .f(Xj).f"(Xt) j kx1) yi 2f(Xl)2 ► rychlost konvergence 3 ► za cenu výpočtu f"\ ► proč tedy nespočítat místo toho další krok N. metody (rychlost konvergence 4)? ► v mnoha případech lze recyklovat části výpočtu, /" je téměř zadarmo ► může i škodit ► pragmaticky ořízneme člen 1 - ^^rlx^ ^° mezí [0.8,1-2] PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 19/52 Seminewtonovské metody ► výpočet derivace nemusí být možný nebo žádoucí ► aproximace f(x + 5)-f(x) není vhodná ► potřebuji 2 x vyhodnocení /, tj. rychost konvergence jen a/2 ► malé 5 - numerická nestabilita ► velké 5 - nepřesnost ► seminewtonovské metody - aproximace derivace „uvnitř" ► metoda sečen ► metoda regula falši PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol < □ ► < [fP ► < -E ► < -E ► -E >0 <\0 20/52 Seminewtonovské metody Metoda sečen ► derivace je aproximována směrnicí spojnice dvou odhadů ► vždy počítá s posledními dvěma odhady ► konverguje obecně rychleji ► zlatý řez (1.618...) ► může porušit separaci kořene PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 21/52 Seminewtonovské metody Regula falši ► důsledně zachovává separaci kořene ► v případě potřeby použije starší odhad PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 22/52 Seminewtonovské metody Regula falsi PA081: Programování numerických výpočtů A. Křenek ► důsledně zachovává separaci kořene ► v případě potřeby použije starší odhad m 2 ^^^^^ / * r * ŕ / t * r * t 4/? A' \ ŕ t t f t 1 * * * w * * í Ě ' * m * * i * * /// i ' * f ' * mi* r 7 p o a * * *■ s * * \* * ' .íí* ► v patologických případech pomalá konvergence Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 22/52 Riddersova metoda ► patologické chování předchozích metod ► jedním z důvodů je proložení přímkou ► idea Riddersovy metody - proložit exponenciálou p(x) = a + be cx ► potřebné 3 body Xo,Xi,X2, omezené na X\ - X() = X'z - X\ = d ► p (x) = 0 v bodě = f(xp) -f(xi) , Mb f(xi) -f(x2) x4=x0 + d-— kde lna . f(x0)-f(x1) b = f(x0) - af(xi) Riddersova metoda ► vyžaduje výpočet dvou logaritmů, příliš náročné ► nahrazeno aproximací x 3 = xq + d u(3 + u2) v(3 + v2) kde u = v = b-l b + l a - 1 a+l ► Ridders, C. A new algorithm for computing a single root of a real continuous function. IEEE Transactions on Circuits and Systems 26: 979-980, 1979. ► X4 vždy spadne do [xo,X2] ► vybereme nejbližší z xo,xi,X2 a dopočítáme třetí do stejné vzdálenosti ► rychlost konvergence V2 ► po krocích kvadraticky, ale vyžaduje dvojí vyhodnocení PA08i: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol Brentova metoda ► půlení intervalu jako bezpečný základ ► proložení inverzní kvadratickou funkcí pro urychlení konvergence ► x jako kvadratická funkce y ► opět tři aktuální body odhadu xi, x2, x$3 výpočet vede na X3 = Xi + Q kde P,Q jsou vyjádřeny z xllx2lx3 a f(x1)1f(x2)1f(x3) elementárními aritmetickými operacemi ► algoritmus hlídá zejména |Q| » 0, jinak se vrací k půlení intervalů ► prakticky zřejmě nejuniverzálnější metoda ► nejsou-li k dispozici derivace ► neřešíme-li speciální případ, kde jiná metoda funguje lépe a/nebo rychleji PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 25/52 Multiplikativní kalkulus ► geometrická derivace fHx) = l\m[f(x + h)Xh ► pro „rozumné" funkce fix) ► platí řada analogií „normálního" kalkulu, zejména Taylorova věta f(x + h)= /(*)/* ix)hr*{x)^ ... /*(n)(x)^/*(n+1)(x + 60^ PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 26/52 Multiplikativní kalkulus ► Newtonova metoda: pro f(x) = l,f*(x) * 1 bude posloupnost Inf(Xj) ln/*(Xi) konvergovat k x ► speciálně pro funkce tvaru egM odstraní exp. složku ► v jistých případech rychlejší praktická konvergence ► kombinace exponenciálních a goniometrických funkcí ► viz http://dx.doi.org/10-1155/2016/8174610 Kořeny polynomů ► speciální případ nelineární rovnice ► lze aplikovat zjednodušená (rychlejší, přesnější) řešení ► silnější sklony ke špatnému chování PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 28/52 Kořeny polynomů ► speciální případ nelineární rovnice ► lze aplikovat zjednodušená (rychlejší, přesnější) řešení ► silnější sklony ke špatnému chování ► špatná podmíněnost, např. Wilkinsonův polynom n i=l ► reálný polynom stupně n má n kořenů Kořeny polynomů Potenciální problémy PA081: Programování numerických výpočtů A. Křenek ► násobné kořeny ► derivace jsou nulové, selhávají Newtonova a seminewtonovské metody ► sudě násobné kořeny nelze separovat ► kořeny mohou být komplexní, co s nimi? Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol Kořeny polynomů Faktorizace kořenů ► pro reálný kořen Pn(x) = (X -Xn)Q_n-i(x) ► pro dvojici komplexních kořenů Pn(x) = (x - (a + íb))(x - (a - íb))Qn_2(x) = (x2 - 2ax + a2 + b2)Qn-2(x) PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 30/52 Kořeny polynomů Faktorizace kořenů ► pro reálný kořen Pn(x) = (X - Xn)Qn-l(x) ► pro dvojici komplexních kořenů Pn(x) = (x - (a + ib))(x - (a - íb))Q_n-2{x) = (x2 - 2ax + a2 + b2)Q_n-2(x) ► známe-li xn, dokážeme spočítat koeficienty Q(x) (qo + (\\x + - xn) = -xnqo + (qo ~ xnqi)x + (qi - xnq2) + ■ ■ ■ a tedy qo = Po x (li = Pi ~ q.i-i n X n PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol ► analogicky opačným směrem, počínaje qn 30/52 Kořeny polynomů Faktorizace kořenů ► dvojitá rekurentní formule, hrozí numerická nestabilita ► vypočteme nepřesné koeficienty Qn-i, použijeme k výpočtu Qn-2, ■■■ ► stabilní chování ► kořen s největší absolutní hodnotou, počítáme od qo ► kořen s nejmenší absolutní hodnotou, počítáme od qn ► „leštění kořenů" ► postupně nalezené kořeny chápeme jako aproximaci ► použijeme v původním P(x) např. do Newtonovy metody ► příliš velká chyba může svést leštění k jinému kořenu (lze ohlídat) PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 31/52 Kořeny polynomů Přímočará metoda ► odchytíme reálný kořen dříve popsanými metodami ► separace metodou pokusu a omylu ► řešení zpravidla Newtonovou metodou ► výpočet derivace polynomu je triviální ► provedeme faktorizaci, opakujeme pro další reálný kořen ► zpravidla včetně „leštění" ► následně faktorizace kvadratických polynomů ► Mullerova metoda - zobecnění metody sečen, aproximace parabolou ► Bairstowova metoda - vede na Newtonovu metodu ve dvou dimenzích PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 32/52 Kořeny polynomů Další metody Laguerre ► ► ► ► ► iterační hledání jednoho kořene včetně komplexních triková manipulace s ln \P(x) | a jeho derivacemi funguje i pro komplexní koeficienty faktorizace a opakované použití iterační krok lze použít k leštění Jenkins-Traub ► propracovaná metoda, základ obecných knihoven ► odvození vyžaduje 4 kapitoly v knize ... Lehmer-Shur ► generalizace separace kořenů na kruhy v komplexní rovině PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 33/52 Polynomy - vyhodnocení výpočetní náročnost PA081: Programování numerických výpočtů A. Křenek P(x) = po + piX + p2X2 + ■ ■ ■ + pnx n ► n násobení xlx ► n násobení ptxi ► n sčítání Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 34/52 Polynomy - vyhodnocení ► výpočetní náročnost P(x) = po + pix + V2X1 + ■ ■ ■ + pnxn ► n násobení xlx ► n násobení pix1 ► n sčítání ► Hornerův algoritmus p[n]=a[n] for (i=n-l;i>=0;i-) p[i]=xp[i+l]+a[i] return p[0] ► pouze 2n operací, výhradně madd ► lze ukázat numerickou stabilitu ► lze přímočaře rozšířit na výpočet derivací PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 34/52 Kořeny polynomů Vlastní hodnoty matic ► vlastní hodnoty matice A jsou kořeny charakteristického polynomu P(x) = det\A-xI\ ► lze zkonstuovat matici / - A = Pm-l Pm 1 0 V m-2 Pm 0 1 Po Pm 0 o V o o o J jejíž charakteristický polynom je právě P(x) = X Ví*1 ► lze aplikovat metody hledání vlastních hodnot ► zpravidla pomalejší ale celkově robustnější PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 35/52 Cirkusové číslo Artisté připravují nové vystoupení. Na začátku se jeden z nich zachytí rukama i nohama v kruhové konstrukci o průměru 180 cm a další konstrukci vykoulí na plošinu. Je třeba, aby se na konci tohoto pohybu původně nejvyšší bod kruhové konstrukce dotýkal plošiny a byl ve stejné výšce jako na začátku. Jaká je potřebná délka a sklon plošiny? PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 36/52 Cirkusové číslo / \ PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 37/52 Cirkusové číslo ► označíme ► R - poloměr kruhu ► a - úhel otočení kruhu (v radiánech) ► j8 - úhel naklonění plošiny ► délka pohybu kruhu po plošině je R a ► z naznačeného čtyřúhelníku odečteme ► a + j8 = tt (další dva úhly jsou pravé) ► tanf = £ ► z toho vyplývá řešená rovnice č 1 tan - =-- 2 TT — p tedy f (p) = (7r-/J)tan| -1 = 0 PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 38/52 Cirkusové číslo separace kořene ► j8 je ostrý úhel, tj. j6 e (0, f) ► řešení je právě jedno ► z geometrické podstaty problému ► formulací rovnic jsme nepřidali falešný kořen ► dáme si práci s exaktním důkazem nebo to rovnou zkusíme ► pro jistotu overíme /(O) = -1 /(7t/2) = 0.5707963 ► numericky nebezpečná by byla až oblast j8 — rr ► vedla by na součin typu 0 ■ oo ► pohybujeme se v bezpečné vzdálenosti ► budeme předstírat, že neumíme spočítat derivaci ► ukážeme metodu sečen, Riddersovu, a Brentovu PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol Cirkusové číslo Metoda puleni intervalu, požadovaná absolutní přesnost 10 n beta f(beta) 1 +0. 0000000 -1. 0000000 2 +1. 5707963 +0. 5707963 3 +0. 7853982 -0. 0240323 4 +1. 1780972 +0. 3119657 5 +0. 9817477 +0. 1544612 6 +0. 8835729 +0. 0679638 7 +0. 8344855 +0. 0226702 8 +0. 8099419 -0. 0005027 9 +0. 8222137 +0. 0111281 15 16 +0.8105171 +0.8104212 +0.0000445 -0.0000467 ► počínaje třetím krokem výpočet vždy ve středu intervalu ► pomalá konvergence Cirkusové číslo Riddersova metoda, požadovaná absolutní přesnost 10 n beta f(beta) 1 +0. 0000000 -1. 0000000 2 +1. 5707963 +0. 5707963 3 +0. 7853982 -0. 0240323 4 +0. 8103685 -0. 0000968 5 +1. 1905824 +0. 3213144 6 +0. 8104702 -0. 0000001 7 +1. 0005263 +0. 1704016 (8) +0. 8104702 4. první iteracní výpočet 5. půlení intervalu mezi 2. a 4. 6. další iterace 7. půlení mezi 5. a 6. 8. výsledek, už je v toleranci Cirkusové číslo Brentova metoda, požadovaná absolutní přesnost 10 -4 n beta f(beta) 1 +0. 0000000 -1. 0000000 2 +1. 5707963 +0. 5707963 3 +1. 0000000 +0. 1699574 4 +0. 7931377 -0. 0165737 5 +0. 8115179 +0. 0009960 6 +0. 8104759 +0. 0000054 7 +0. 8104259 -0. 0000422 (8) +0. 8104759 ► počítá primárním způsobem, nedošlo na bezpečné půlení intervalů ► srovnatelně rychlá konvergence s Riddersem ► při větší požadované přesnosti se liší ± o jeden krok PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 42/52 Cirkusové číslo ► všechny metody se v rámci požadované tolerance shodují ► rychlost konvergence dle očekávání ► zdvojnásobení požadované přesnosti zdvojnásobí počet kroků půlení intervalů ► Riddersova metoda naroste ze 7 na 9 ► Brentova dokonce jen ze 7 na 8 ► lepší výsledky než teoretická rychlost konvergence \/2 Cirkusové číslo ► délka pohybu po plošině je Ra = R(tt-P) = 2.098 ► zbývá spočítat část plošiny od země k bodu dotyku 7T — j8 R tan—-— = — tedy i = R tan((7r-0)/2) = 0.3861 ► nehrozí významné numerické nepřesnosti ► celková délka plošiny je 2.484 m, sklon 46.44' PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 44/52 Obchodní případ Americká garážová firma zahajuje montáž počítačů. Kolik jich musí v prvním roce prodat, aby byla zisková? PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 45/52 Obchodní případ Americká garážová firma zahajuje montáž počítačů. Kolik jich musí v prvním roce prodat, aby byla zisková? ► uvažované náklady ► fixní jednorázové (pořízení vybavení): $20000 ► fixní roční (nájem, web, ...): $15000 ► variabilní (komponenty, hodinová sazba práce, ...): $625n ► semivariabilní (náročnější logistika atd.) $30n15 ► předpokládané příjmy ► prodej počítačů: $1500n ► slevy (množstevní, VIP zákazníci, ...): -$10n15 PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 45/52 Obchodní případ ► řešená rovnice f (n) = TC(n) - TS(n) = 35000 - 875n + 40n15 = 0 ► první dva členy - lineární funkce ► prochází bodem (0,35000) ► osu x protíná v 35000/875 = 40 ► člen 40n15 způsobí „prohnutí" nahoru ► posune kořen z bodu 40 dál ► přidá druhý kořen pro vyšší n ► zisku dosahujeme v oblasti mezi těmito kořeny PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 46/52 Obchodní případ ► separace 1. kořene ► podle předchozí úvahy /(40) > 0, skutečná hodnota 10119 ► zkusíme /(80) = -6378, vyhovuje ► separace 2. kořene ► předchozí /(80) = -6378 ► hrubý odhad zanedbáním absolutního člene -875n + 40n15 = 0 tedy ► zkusíme /(500) = 44713, vyhovuje ► numerická stabilita v dané separaci kořenů ► součty/rozdíly čísel s minimálním řádovým rozdílem ► reálně nás zajímá přesnost na jednotky ► žádný potenciální problém PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol Obchodní případ ► vypočtené kořeny 62.7 a 384.0 ► výsledky v požadované přesnosti shodné všemi metodami ► chování jednotlivých metod (počty iterací) přesnost 10 -i 10 -4 kořen 1 2 1 2 půlení 11 15 21 25 Ridders 5 9 7 11 Brent 5 9 7 10 ► umělá přesnost 10"4 - zdvojnásobení počtu platných číslic ► dokládá očekávanou rychlost konvergence ► lineární pro půlení intervalů "2 pro ostatní metody PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 48/52 Shrnutí ► řešíme nelineární rovnice tvaru F(x) = G(x) iteračními metodami ► redukce na hledání kořenů f(x) = 0 ► separace kořenů je klíčová ► nalezení intervalu, ve kterém leží právě jeden ► metoda půlení intervalů ► robustní ale pomalá ► Newtonova metoda ► rychlá konvergence ► vyžaduje derivaci a dobrý počáteční odhad ► seminewtonovské metody ► vnitřní aproximace derivace ► metoda sečen, regula falši ► pokročilé metody ► aproximace funkce složitější křivkou ► Ridders, Brent ► speciální metody pro polynomy PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 49/52 Domácí úkol František poseděl s přáteli v restauraci a vrací se poněkud nejistým krokem domů. Cesta měří 2400 m a vede poblíž meandrujícího potoka. Spočítejte, zda a jak daleko od začátku cesty František do potoka spadne. PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 50/52 Domácí úkol František poseděl s přáteli v restauraci a vrací se poněkud nejistým krokem domů. Cesta měří 2400 m a vede poblíž meandrujícího potoka. Spočítejte, zda a jak daleko od začátku cesty František do potoka spadne. K dispozici máte: ► double frantisek(double x, double *df) ► pro 0 < x < 2400 na vrací Františkovu souřadnici y ► není-li df == NULL, spočítá také derivaci dy/dx ► při každém běhu programu počítá jinak! ► double potok(double x, double *dp) ► totéž pro potok, počítá vždy stejně PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 50/52 Domácí úkol 20 15 10 5 0 -5 -10 -15 -20 0 50 100 150 200 PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol 51/52 Domácí úkol ► Nejsou přístupné zdrojové kódy ► nevíte tedy přesně, jak se funkce chovají ► K separaci kořene využijte „přirozené" charakteristiky ► meandry potoka jsou dlouhé v řádu desítek metrů ► alkoholem ovlivněný krok se může vychýlit do strany i každý metr ► obě funkce jsou spojité ► experimentujte s hustotou vzorkování ► Počítejte s uměle vysokou přesností 10~12 v y ► Implementujte metody bez derivace i s derivací, srovnejte náročnost ► funkce long cena() vrací počet vyhodnocení obou funkcí, derivace se počítá 2 x ► Termín odevzdání 31.3., 1 bod ke zkoušce PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewtonovské metody Sečny Regula falši Riddersova metoda Brentova metoda Polynomy Příklad Shrnutí Domácí úkol