PA081: Programování numerických ^ - o vypočtu 2. Nelineární rovnice o jedné neznámé Aleš Křenek jaro 2016 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 1/54 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í 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí Metoda prosté iterace ► rovnice ve tvaru x = f(x) ► řešením je pevný bod funkce / ► počítáme opakovaně Xí+i = f(xt) ► až k dosažení kritéria konvergence PA081: Programovaní 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 4/54 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí Metoda prosté iterace Pro železniční příklad PA081: Programování numerických výpočtů A. Křenek ► je zobrazení a »- *f síikx kontrakce? l/L ► postačující podmínka a V x e K: f'(x) < 1 tedy cos a < — a d ► | siná = a bude mít řešení v (arccos f, 7) Formulace problému a příklad Prostá iterace Separace kořenů Půlení intervalů Newton Seminewto no vs ké metody Sečny Regula falši Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 5/54 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 ► ^ sinex = 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) 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(CLi)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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 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 / h 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 8/54 Separace kořenů Nespojitá funkce ► 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ř. -±rc ► některé metody najdou „kořen" v c ► snadno identifikovatelný problém 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 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 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 ► 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí < □ ► < rS1 ► < ► < ► -E f)0>O 11/54 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 12/54 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 13/54 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí Newtonova metoda Konvergence ► označíme x* kořen, a e i odchylky X{ - x x" + ei+1 = Xi+i = x tedy ei+i = et -Taylorův rozvoj f(Xj) 1 f'(Xi) = X* + €i - f(Xj) 1 f'(Xi) 0 = f(X*)=f(Xi)-€if(Xi) + 2! kde g [0, cí] ► podělením f'(xt) a dosazením dostaneme ei+l = ~ei 2/"(&) 2f'(Xi) 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí ► kvadratická konvergence ► v každém kroku se zdvojnásobí počet platných číslic 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 * w + m 1 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí = ► 16/54 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 17/54 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 \/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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 18/54 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 /(v) ✓ # + m * Ě + m * § * g i > i * * v t 1 ' M V-/ M * 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí Seminewtonovské metody Regula falši ► důsledně zachovává separaci kořene ► v případě potřeby použije starší odhad m 2 / * r * ŕ * f f A' \ * 1 * * * í • # * * m * * i * * /// i ' * i * * r 7 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 20/54 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 20/54 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) x4 = x0 + d-— kde lna . f(x0)-f(x1) b = f(x0) - af(xi) Riddersova metoda ► vyžaduje výpočet dvou logaritmu, 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í 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 22/54 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 xi,x2,x3 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 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 ► násobné kořeny ► derivace jsou nulové, selhávají Newtonova seminewtonovské metody ► sudě násobné kořeny nelze separovat ► kořeny mohou být komplexní, co s nimi? < □ ► < rgi ► 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 26/54 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 - ib))Q_n-2{x) = (x2 - 2ax + a2 + b2)Qn-2(x) ► známe-li xn, dokážeme spočítat koeficienty Q(x) (qo + 4ix + ...)(x - xn) = -xnqo + (qo - xnqi)x + (q± - xnq2) + ■ ■ ■ a tedy qo = Po x (li = Pí - qt-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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí analogicky opačným směrem, počínaje qn 26/54 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 27/54 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 28/54 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí Polynomy - vyhodnocení ► výpočetní náročnost P(x) = po + P\X + p2x2 + ■ ■ ■ + pnx ► n násobení xlx ► n násobení ptx1 ► n sčítá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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí Polynomy - vyhodnocení ► výpočetní náročnost P(x) = po + pix + P2X2 + ■ ■ ■ + pnxn ► n násobení xlx ► n násobení pix1 ► n s cit ani ► Hornerův algoritmus p[n]=a[n] foř (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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 30/54 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 31/54 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí Montáž ložiska K sestavení kluzného ložiska uložení ocelové mostní konstrukce je třeba na místě zasunout čep o vnějším průměru 313 mm do náboje o stejném vnitřním průměru. Sestavení probíhá s podchlazeným čepem, díky tepelné roztažnosti materiálu je menší a ložisko lze sestavit. Bezpečné zmenšení průměru je 0.37 mm. Na jakou teplotu je třeba čep podchladit, předpokládáme-li teplotu prostředí 20°C? 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 33/54 Montáž ložiska ► koeficient tepelné roztažnosti oceli a závisí na teplotě ► výpočet tedy není zcela triviální ► v rozmezí ±200°C je poměrně přesně vyjádřen funkcí a(T) =aT2 + bT + C kde a = -8.27 x ÍO-11,*? = 2.21 x l(r8,c = 1.17 x 1(T5 ► koeficienty získány regresí (proložením křivky) experimentálních dat ► metodu regrese budeme probírat později Montáž ložiska ► poloměr po zchlazení počítáme pro AT — 0 rx = r0(l - ATa(T0))(l - ATa(T0 - AT))... = r0(l - ATa(To) - ATa(T0 - AT) + (AT)2a(T0)(x(To-AT))... ► vedlo by na netriviální diferenciální rovnici ► protože a(T) <*c 1, AT < T2 ► nemají fyzikální význam ► T+ - čep se při vysokých teplotách začne smršťovat ► T- - při nízkých teplotách se začne znovu roztahovat ► obojí důsledkem aplikace regresní křivky mimo rozsah hodnot, ze kterých byla určena ► skutečný kořen v intervalu [Ti, T2] ► funkce je spojitá, kořen je právě jeden ► dokonalá separace ► umíme jednoduše počítat derivaci ► navíc je na (Ti, t2) derivace nenulová ► Newtonova metoda je ideální Montáž ložiska ► numerická stabilita výpočtu f (T) a f (T) ► rozsah T řádově ve stovkách ► členy polynomu v řádech 10~5-10~2 ► raději použijeme typ doubl e ► potenciálni problém Newtonovy metody f (T) ► nenastane díky známým vlastnostem funkce Montáž ložiska ► Průběh výpočtu funkcí rtnewt ► požadovaná absolutní přesnost 10~3 ► nemá smysl více, vstupní data (koeficienty a, b, c) jsou méně přesné ► a jak budeme při stavbě mostu chladit čep na takhle přesnou teplotu? ► v následujícím kroku už bylo dosaženo požadované přesnosti ► výsledek je -90.056°C 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? 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 42/54 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 Ra ► 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 43/54 Cirkusové číslo separace kořene ► P je ostrý úhel, tj. j8 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 44/54 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 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 a/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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 49/54 Obchodní případ Americká garážová firma zahajuje montáž počítačů. Kolik jich musí v prvním roce prodat, aby byla zisková? 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, . 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 ► ► ► ► ): $625n ► ► 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 50/54 Obchodní případ ► řešená rovnice f (n) = TC(n) - TS(n) = 35000 - 875n + 40n1 ► 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 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 53/54 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 Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 54/54