PA081: Programování numerických výpočtů
3. Nelineární rovnice o jedné neznámé
Aleš Křenek
jaro 2011
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í
► 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
PA081: Programování numerických
výpočtů
A. Křenek
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 = R(l - cos a)
dosazením a jednoduchou úpravou
d .
— sin a = a a
hledáme pevný bod funkce
kandidát na metodu prosté iterace
3/34
Metoda prosté iterace
rovnice ve tvaru x = f (x) řešením je pevný bod funkce / ► počítáme opakovaně xi+i = f(xt) ► 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ě xi+i = f(xt)
► až k dosažení kritéria konvergence
► kdy a proč to funguje? - věta o pevném bodě
Je-li pro ře| funkce /: K ->■ K kontrakce, tj. existuje q e (0,1) tak, že Mx,y <eK: \ f(x) -f(y)\ < q\x - y\, potom existuje x* e K takové, že pro libovolné xo e K x* limitou posloupnosti xt+i = f(xt).
► idea důkazu
|xi+i-x*| = \f(xi)-f(x*)\< \Xi-X*\
Metoda prosté iterace
Pro železniční příklad
► je zobrazení a    ^ sin a kontrakce?
► postačující podmínka
Vx e K: /'(%) < 1   tedy   cos a < ^
a
► f sin a = a bude mít řešení v (arccos j, y)
PA081: Programování numerických
výpočtů
A. Křenek
5/34
Metoda prosté iterace
Pro železniční příklad
► je zobrazení a    ^ sin a kontrakce?
► postačující podmínka
Vx e K: /'(%) < 1   tedy   cos a <
a
► f sin a = a bude mít řešení v (arccos j, 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 [ai, bi]
► v každém [ai, bi\ má funkce právě jeden kořen
► f(cu)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 [at,bi]
► stanovíme konkrétní podmínku ukončení
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
PA081: Programování numerických
výpočtů
A. Křenek
7/34
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
8/34
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á
- některé metody najdou „kořen" v c
► snadno identifikovatelný problém
Separace kořenů
Patologické případy
	U u
\              \ \ e         f c	•     1              1 \ X\ d         a X2  X3 b
PA081:
Programování numerických výpočtů
A. Křenek
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 j atd., je problém zpravidla někde jinde
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 j 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
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 j 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 -* ±<x> opačné zneménko
11/34
Separace kořenů
Algoritmus „look outward"
PA081: Programování numerických
výpočtů
A. Křenek
int zbrac(float (*func)(float), float *xl, float *x2) {
i nt j ;
float fl,f2; fl=(*func)(*xl); f2=(*func)(*x2); for (j=l;j<=NTRY;{
if (fl*f2 < 0.0) return 1; if (fabs(fl) < fabs(f2)) fl=(*func)(*xl += FACT0R*(*xl-
rx2));
el se
f2=(*func)(*x2
}
return 0;
+= FACT0R*(*x2-*xl));
12/34
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
/(^)/(a)<0 f(^)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
Konvergence
► je-li požadovaná přesnost e, je třeba
Metoda půlení intervalů
PA081: Programování numerických
výpočtů
A. Křenek
n = log2
b - a
e
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ů man tisy
► 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
14/34
Metoda půlení intervalů
Konvergence
► je-li požadovaná přesnost e, je třeba
b - a n = log2 ——
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ů man tisy
► 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á
15/34
Newtonova metoda
► také Newton-Raphsonova metoda
► odvozena z Taylorova rozvoje x*
f(xi + 5i) =f(Xi)+f(x)5
► zanedbáme vyšší derivace
e     . f{Xi)
Oi -
fix i)
► opakujeme iterační krok
f(Xi)
Xi+l — Xi
f'(Xi)
► nutno spočítat i derivaci
Newtonova metoda
Konvergence
► označíme x* kořen, a ei odchylky xi - x*
* f{Xi) o f(Xi)
x* + ei+i = xf+i = xt- 4tV^T = x* + et
tedy ei+i = eť -
► Taylorův rozvoj
0 = /(**) = f(xi) - ej'ixi) + e2if"^i]
kde ŠíG [0,eť]
► podělením/'(xí) a dosazením dostaneme
£í+l ~ ~£í n r,
f" {li)
2f'{xi)
► 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
PA081: Programování numerických
výpočtů
A. Křenek
		i	
			
prakticky nepoužitelná sama o sobě ► nutno kombinovat s jinými metodami vhodná k „vyleštění" nahrubo nalezeného kořene
Newtonova metoda
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<=!MAX;{
(*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říliš mnoho iterací */ return 0;
19/34
Seminewtonovské metody
PA081: Programování numerických
výpočtů
A. Křenek
► výpočet derivace nemusí být možný nebo žádoucí
► aproximace
není vhodná
► potřebuji 2 x vyhodnocení /, tj. rychost konvergence jen -J2
► malé 5 - numerická nestabilita
► velké 5 - nepřesnost
► seminewtonovské metody - aproximace derivace „uvnitř"
► metoda sečen
► metoda regula falši
fix)
fix+ 5)-f{x)
ô
20/34
Seminewtonovské metody
Metoda sečen
► derivace je aproximována směrnicí spojnice dvou odhadů
PA081: Programování numerických
výpočtů
A. Křenek
vždy počítá s posledními dvěma odhady konverguje obecně rychleji ► zlatý řez (1.618...)
může porušit separaci kořene
21/34
Seminewtonovské metody
Regula falši
► důsledně zachovává separaci kořene
► v případě potřeby použije starší odhad
Seminewtonovské metody
Regula falsi
► v patologických případech pomalá konvergence
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
potřebné 3 body xq, x\,xz, omezené na
x\ - Xq = X2 - x\ = d
► p (x) = 0 v bodě
p(x) = a + be'
f(x0) - f{xi) f(xi)-f(x2) f(x0) - f(Xl) f(x0) - af(xi)
X4 = Xq + d
\wb lna
a
kde
Riddersova metoda
► vyžaduje výpočet dvou logaritmů, příliš náročné
► nahrazeno aproximací
%3 = Xq + d
u(3 + u2) v(3 + v2)
kde
u
v
b
b + l a - 1
a + 1
► 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.
%4 vždy spadne do [xo,X2]
► vybereme nejbližší z xo, x\,X2 a dopočítáme třetí do stejné vzdálenosti
► rychlost konvergence \Ti
► po krocích kvadraticky, ale vyžaduje dvojí vyhodnocení
24/34
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 x\,X2,x^, výpočet vede na
P
x-i = xi + —
kde P, Q jsou vyjádřeny zx1,x2,x3af(x1),f(x2),f(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
Domácí úkol
► implementujte řešení železničního příkladu probranými základními metodami
► prostá iterace, Newton, půlení intervalů, sečny a/nebo regula falši
► vstupem je i požadovaná přesnost výsledku (vzdálenost středu od spojnice)
► testujte na více různých vstupech
► pro každou metodu nalezněte co nejlepší separaci kořene/počáteční odhad
► bezpečně - bude konvergovat
► co nejlepší konvergence
► odevzdejte
► implementaci ve svém oblíbeném programovacím jazyce
► zhodnocení konvergence metod pro tento příklad
► přiměřené komentáře k volbě počátečních hodnot
► hodnocení: 2 body
► akceptuji týmové řešení (skupinky 2-3)
Domácí úkol
► numericky stabiní řešení i pro velmi málo odlišná a a d v zadání
► nápověda: důvěryhodnost metody se projeví srovnáním vypočteného poloměru zakřivení
► hodnocení: 2 body
Domácí úkol
► numericky stabiní řešení i pro velmi málo odlišná a a d v zadání
► nápověda: důvěryhodnost metody se projeví srovnáním vypočteného poloměru zakřivení
► hodnocení: 2 body
► připravte krátkou prezentaci nej zajímavějších fenoménů
► hodnocení: f bod
► termín: 31.3.2011
Kořeny polynomů
PA081: Programování numerických
výpočtů
A. Křenek
► 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í
28/34
Kořeny polynomů
A. Křenek
► 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
Wn(x) = Y\(X-Í) í=l
► reálný polynom stupně n má n kořenů
PA081: Programování numerických výpočtů
28/34
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?
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))Q_n-2M
= (x2 - 2ax + a2 + b2)Q_n-2(x)
PA081: Programování numerických
výpočtů
A. Křenek
30/34
Kořeny polynomů
Faktorizace kořenů
► pro reálný kořen
Pn(x) = (X - Xn)Q.n-l(x)
► pro dvojici komplexních kořenů
Pn(x) = (x - (a + íb))(x - (a - íb))Q_n-2(x)
= (x2 - 2ax + a2 + b2)Qn
► známe-li xn, dokážeme spočítat koeficienty Q(x)
(q0 + qix + ... )(x - xn)
= ~xnqo + (qo ~ xnqi)x + (q1 - xnq2'.
a tedy
n _   Po      n _ Pí ~ (ji-i
qo---qi--
Xn Xn
► analogicky opačným směrem, počínaje qn
Kořeny polynomů
Faktorizace kořenů
► dvojitá rekurentní formule, hrozí numerická nestabilita
► vypočteme nepřesné koeficienty Q_n-\, použijeme
k Výpočtu Qm-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)
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
32/34
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
33/34
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-1	Pm-2	Po
Pm	Pm	Pm
1	0	0
0	1	0
0	0	0
jejíž charakteristický polynom je právě P(x) = X ptx1 ► lze aplikovat metody hledání vlastních hodnot
► zpravidla pomalejší ale celkově robustnější