PA081: Programování numerických výpočtů 3. Nelineární rovnice o jedné neznámé Aleš Křenek jaro 2013 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 Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Re gula falši Domácí úkol Polynomy Shrnutí □ s 4 : 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). Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Re gul a falši Domácí úkol Polynomy Shrnuti x i+l X \f(Xi)-f(x*)\ < \Xi-X* Metoda prosté iterace Pro železniční příklad ie zobrazení a >— - sincx kontrakce? postačující podmínka a Vx g K: f'(x) < 1 tedy cos a < — a srna = a bude mít řešení v (arccos |, ^) Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Re gul a falši Domácí úkol Polynomy 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 Re gul a falši Domácí úkol Polynomy Shrnutí □ 4 ► 4 š ► < -E ► 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 Re gul a falši Domácí úkol Polynomy Shrnutí Separace kořenů Nespojitá funkce PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad není-li / spojitá, ale je ohraničená ► „kříží" osu x v bodě nespojitosti ► numericky nerozlišitelné od kořene Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Re gul a falši Domácí úkol Polynomy Shrnutí Separace kořenů Nespojitá funkce PA081: Programování numerických výpočtů A. Křenek Formulace problému a příklad není-li / spojitá, ale je ohraničená ► „kříží" osu x v bodě nespojitosti ► numericky nerozlišitelné od kořene není-li ani ohraničená K naPř- ^ ► některé metody najdou „kořen" v c *■ snadno identifikovatelný problém Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Re gul a falši Domácí úkol Polynomy Shrnutí □ 4 ► 4 ■= * < ► 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 Re gul a falši Domácí úkol Polynomy Shrnuti 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 Re gul a falši Domácí úkol Polynomy 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" (NRC) ► 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 Re gul a falši Domácí úkol Polynomy 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" (NRC) ► 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" (NRC) ► 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 Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Re gul a falši Domácí úkol Polynomy Shrnutí □ 4 : Separace kořenů Algoritmus „look inward" void zbrak(float (*fx) (float), float xl, float x2, int n, float xbl[], float xb2[], int *nb) { int nbb,i; float x,fp,fc,dx; nbb=0; dx=(x2-xl)/n; fp=(*fx)(x=xl); for (i=l; i<=n; i++) { fc=(*fx)(x += dx); if (fc*fp <= 0.0) { xbl[++nbb]=x-dx; xb2[nbb]=x; if(*nb == nbb) return; } fp=fc; } *nb = nbb; Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Re gula falši Domácí úkol Polynomy Shrnutí Separace kořenů Algoritmus „look outward" int { zbrac(float (*func) (float) , float *xl, float *x2) int j; float fl,f2; fl=(*func)(*xl); f2=(*func)(*x2); for (j =1;j<=NTRY;j ++) { if (fl*f2 < 0.0) return 1; if (fabs(fl) < fabs(f2)) fl=(*func)(*xl += FACTOR*(*xl-*x2)); else f2=(*func)(*x2 += FACT0R*(*x2-*xl)); } return 0; 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 2 )J J \ 2 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 Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Re gul a falši Domácí úkol Polynomy Shrnutí □ 4 ► 4 š ► < -E ► Metoda půlení intervalů Konvergence ► je-li požadovaná přesnost e, je třeba n log; b - a iteračních kroků 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\ Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Re gul a falši Domácí úkol Polynomy Shrnuti kde e je přesnost daného datového typu metoda je robustní ale relativně pomalá Metoda půlení intervalů float rtbis(float (*func) (float), float xl, float x2, float xacc,int *iter) { i nt j ; float dx,f,fmid,xmid, rtb; f=(*func)(xl) ; fmid=(*func) (x2) ; rtb = f < 0.0 ? (dx=x2-xl,xl) : (dx=xl-x2,x2); for (j=l;j<=JMAX;j++) { fmid=(*func)(xmid=rtb+(dx *= 0.5)); if (fmid <= 0.0) rtb=xmid; if (fabs(dx) < xacc || fmid ==0.0) { '•iter = j; return rtb; } } Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Re gul a falši Domácí úkol Polynomy Shrnuti □ 4 (5 ► 4 ■= * < ► Newtonova metoda ► také Newton-Raphsonova metoda ► odvozena z Taylorova rozvoje x* = xt + 5i f"(Xi)5\ f(xi + 5i) = f{Xi) + f'{Xi)5i + 2! zanedbáme vyšší derivace 5i f(Xj) f(xt) opakujeme iterační krok fJXj) Xi + 1 — Xi „ f (Xi) Formulace problému a příklad Separace kořenů Půlení intervalů Newton Seminewtonovsk metody Re gul a falši Domácí úkol Polynomy Shrnutí nutno spočítat i derivaci Newtonova metoda Konvergence ► označíme x* kořen, a ei odchylky xt - x* x* + ei+1 = xi+1 = Xi- = x* + et fiXi) f(Xi) 1 f'(Xt) f(Xj) f(Xi) Taylorův rozvoj tedyeŕ+1=eŕ-$S) 0 = fix*) = f(Xi) - EifiXi) + ^ (gr) 2! kde &e [O.ei] podělením f'(xi) a dosazením dostaneme .2 f "(li) ei+i kvadratická konvergence 'l2f'(Xi) v každém kroku se zdvojnásobí počet platných číslic □ ►