PA081: Programování numerických výpočtů 2. Nelineární rovnice o jedné neznámé Aleš Křenek jaro 2015 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í Xq ► 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í Xq ► 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 xq ► 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 2ti 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 R sin a = a R(l COS Oř) dosazením a jednoduchou úpravou d sin a a a hledáme pevný bod funkce kandidát na metodu prosté iterace □ s 4 : Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Re gula 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ě xt+i = f {x i) ► až k dosažení kritéria konvergence Metoda prosté iterace rovnice ve tvaru x = f (x) řešením je pevný bod funkce / počítáme opakovaně xt+i = f {x i) ► až k dosažení kritéria konvergence kdy a proč to funguje? - věta o pevném bodě Je-li pro K c K funkce /: K —■ K kontrakce, tj. existuje q g (0,1) tak, že \/x,y g K: \f{x) - f{y)\ < q\x - y\, potom existuje x* g K takové, že pro libovolné xq g K je x* limitou posloupnosti xt+i idea důkazu f(Xi). x ŕ+1 x \f(Xi)-f(x*)\ < \Xi-X* Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Regula falši Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí Metoda prosté iterace Pro železniční příklad je zobrazení a >— - sincx kontrakce? postačující podmínka a V x G K: f (x) < 1 tedy cos cx < — a sincx = a bude mít řešení v (arccos §, y) □ s 4 : Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Regula falši Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí Metoda prosté iterace Pro železniční příklad ► ie zobrazení cx >— 4 sincx kontrakce? *■ postačující podmínka Vx G K: f (x) < 1 tedy cos a < —3 a ► ^ sincx = a bude mít řešení v (arccos §, y) ► 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 [au bt\ *■ v každém [au bi\ má funkce právě jeden kořen ► f(at)f(bt) < 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 [cii,bi] stanovíme konkrétní podmínku ukončení Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Regula falši Polynomy Příklady Montáž ložiska Cirkusové čislo Obchodní pfipad 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 / b 1 m 1 Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody 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 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ř- lh ► některé metody najdou „kořen" v c ► snadno identifikovatelný problém Separace kořenů Patologické případy U u l \ l e f c X\ d a x2 x i b Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody 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 Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody 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 Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody 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 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 4 □ M0M1 M 1 » 1 je Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Regula falši Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 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 /(^)/ n násobení ptx1 *■ n sčítání Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Re gul a falši Polynomy Příklady < □ ► 4 ► 4 = > 4 M > -OQ.O 30/54 Polynomy - vyhodnocení ► výpočetní náročnost P(x) = Po + Pix + P2X2 + ■ ■ ■ + ► n násobení xlx ► n násobení pix1 ► n sčítání ► Hornerův algoritmus ► pouze 2n operací, výhradně madd ► lze ukázat numerickou stabilitu ► lze přímočaře rozšířit na výpočet derivací p[n]=a[n] for (i=n-l;i>=0;i-) p[i]=xp[i+l]+a[i] return p[0] Kořeny polynomů Vlastní hodnoty matic vlastní hodnoty matice A jsou kořeny charakteristického polynomu P(x) = det \A - xl\ lze zkonstuovat matici Pm-1 Pm 1 A 0 V o Pm-2 Pm 0 1 0 Po \ Pm 0 o o / jejíž charakteristický polynom je právě P(x) = X PiX1 lze aplikovat metody hledání vlastních hodnot ► zpravidla pomalejší ale celkově robustnější Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Re gul a falši Polynomy Příklady Montáž ložiska 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 Regula falši Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 32/54 Príklady 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.3 7 mm. Na jakou teplotu je třeba čep podchladit, předpokládáme-li teplotu prostředí 20°C? Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Regula falši Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí Montáž ložiska ► koeficient tepelné roztažnosti oceli cx 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) = aTz + bT + C kdea = -8.27 x ÍO"11,*? = 2.21 x l(T8,c = 1.17 x 10" ► 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 - Ar«(7b))(l - ATa(T0- AT))... = r0(l - ATa(To) - ATa(T0 - AT) + (AT)za(T0)a(T0 - AT))... ► vedlo by na netriviální diferenciální rovnici ► protože a(T) <<; 1, AT <<; 1, lze nelineární členy zanedbat ► dostáváme zjednodušení rx = r0(l - á(Tx)) a budeme řešit rovnici kde á(Tx r To (x(T)dT ro Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Regula falši Polynomy Příklady Montáž ložiska Cirkusové čislo Obchodní případ Shrnutí Montáž ložiska ► řešená funkce f(Tx) = á(Tx -(T3 ro ro b _3T Tl) - potřebujeme ► separaci „správného" kořene ► zvolit vhodnou numerickou metodu Montáž ložiska řešená funkce f(Tx) = á(Tx -(T3 ro ro b a ^ b ^7 3T +2T c(T0 potřebujeme ► separaci „správného" kořene ► zvolit vhodnou numerickou metodu koeficient u T% je kladný (a bylo záporné) pro Tx pro Tx o také/(Tx) --oo také f(Tx) Montáž ložiska ► derivace řešené funkce f(Tx) = -aTx-bTx-c * f(Tx) má jedno lokální minimum a jedno maximum ► v bodech, kde/'(Tx) =0 ► numericky stabilní výpočet? ► první přednáška;-) ► b = -2.21 x 10-8, VB = 6.60 x 10-8 je OK ► při nepřesném určení by mohla Newtonova metoda zabloudit, viz dále ► vypočtené extrémy Ti = -265.544 /(7\) = 0.000867609 T2 = 532.775 f(T2) = -0.00614507 Montáž ložiska ► dva „nesmyslné" kořeny T- < Ti a T+ > 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 [ T\, T2 ] ► funkce je spojitá, kořen je právě jeden ► dokonalá separace ► umíme jednoduše počítat derivaci ► navíc je na (Ti, Ti?) 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 double ► 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? f (T) f (T) 133.615 -66.6454 -89.1 -90.0546 -90.0564 -0.00263873 -0.000221398 -8.66255e-06 -1.68088e-08 -6.3964e-14 -1.31765e-05 -9.85981e-06 -9.07435e-06 -9.03911e-06 -9.03904e-06 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? Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Regula falši Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí Cirkusové číslo Cirkusové číslo označíme ► R - poloměr kruhu ► a - úhel otočení kruhu (v radiánech) ► - úhel naklonění plošiny délka pohybu kruhu po plošině je Ra z naznačeného čtyřúhelníku odečteme ► a + = n (další dva úhly jsou pravé) tan Ě. - 2 Ra z toho vyplýva řešená rovnice tan 1 TT — tedy f (P) = (rr-^)tan^ - 1 = 0 Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Regula falši Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí Cirkusové číslo separace kořene ► P je ostrý úhel, tj. P 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 ověříme /(O) 1 f(n/2) = 0.5707963 numericky nebezpečná by byla až oblast fi — tt ► vedla by na součin typu 0 ■ co ► 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 Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Regula falši Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí Cirkusové číslo Metoda puleni intervalu, požadovaná absolutní přesnost 10~4 beta f(beta) 1 +0. .0000000 -1. .0000000 2 +1. .5707963 +0. .5707963 3 +0. .7853982 -0. .0240323 4 +1. .1780972 +0. .3119657 +0. .9817477 +0. .1544612 6 +0. .8835729 +0. .0679638 +0. .8344855 +0. .0226702 8 +0. .8099419 -0. .0005027 9 +0. .8222137 +0. .0111281 +0. .8105171 +0. .0000445 +0. .8104212 -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 f(beta) -1.0000000 +0.5707963 -0.0240323 -0.0000968 +0.3213144 -0.0000001 +0.1704016 4. první iterační 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 n beta 1 +0.0000000 2 +1.5707963 3 +0.7853982 4 +0.8103685 5 +1.1905824 6 +0.8104702 7 +1.0005263 (8) +0.8104702 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 Separace kořenů Půlení intervalů Newton Seminewtonovsk metody 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 y/2 Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Regula falši Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí Cirkusové číslo délka pohybu po plošině je Ra = R(n - fi) = 2.098 zbývá spočítat část plošiny od země k bodu dotyku n-P R . , R tan ■ 7 tedy l tan((7r-0)/2) nehrozí významné numerické nepřesnosti celková délka plošiny je 2.484 m, sklon 46.44° 0.3861 Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Regula falši Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí Obchodní případ Americká garážová firma zahajuje montáž počítačů. Kolik jich musí v prvním roce prodat, aby byla zisková? Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Regula falši Polynomy Příklady Montáž ložiska Cirkusové čislo Obchodní případ Shrnutí 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, ...): $ f 5000 ► variabilní (komponenty, hodinová sazba práce, ...): $62 ► semivariabilní (náročnější logistika atd.) $30n15 ► předpokládané příjmy ► prodej počítačů: $1500n ► slevy (množstevní, VIP zákazníci, ...): -$1 On1-5 Obchodní případ ► řešená rovnice f (n) = T C (n) - T S (n) = 35000 - 875n + 40ti1 ► první dva členy - lineární funkce ► prochází bodem (0,35000) ► osu x protíná v 35000/875 = 40 ► člen 40n1-5 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 + 40nL5 = 0 tedy n = —- = 478.51 ► 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 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ů ► y/2 pro ostatní metody Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Regula falši Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí 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 Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Regula falši Polynomy Příklady Montáž ložiska Cirkusové číslo Obchodní případ Shrnutí