Hledání minima funkce jedné proměnné

MASARYKOVA UNIVERZITA PŘÍRODOVĚDECKÁ FAKULTA

Hledání minima funkce jedné proměnné

BAKALÁŘSKÁ PRÁCE

Jaroslav Urbánek

Brno, jaro 2006

Prohlášení Prohlašuji, že tato práce je mým původním autorským dílem, které jsem vypracoval samostatně. Všechny zdroje, prameny a literaturu, jež jsem při vypracování používal nebo z nich čerpal, v práci řádně cituji s uvedením úplného odkazu na příslušný zdroj.

V Brně, 23.5.2006 Jaroslav Urbánek

Vedoucí práce: Prof. RNDr. Ivana Horová, CSc.

ii

Poděkování Rád bych na tomto místě poděkoval vedoucímu této práce, prof. RNDr. Ivaně Horové, CSc., za věnovaný čas, rady a připomínky, které mi pomohly k jejímu vypracování.

iii

Shrnutí Práce představuje a rozebírá pět základních numerických metod pro hledání minima funkce jedné proměnné. Vedle hlavní myšlenky prezentuje každou z metod na stejném příkladu pro možnost srovnání. Důraz je kladen na představení metod a jejich použití v praxi. Celkové zhodnocení kladů a záporů jednotlivých metod je uvedeno v sekci Závěr. Součástí práce je CD s programy demonstrujícími jednotlivé metody na konkrétních příkladech.

iv

Klíčová slova minimum, funkce, zlatý řez, kvadratická interpolace, bisekce, Bolzanova metoda, Newtonova metoda, metoda tečen, regula falsi, optimalizace

v

Obsah 1 Úvod

2

2 Pojmy a předpoklady

3

3 Přímé metody 3.1 Zlatý řez . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Kvadratická interpolace . . . . . . . . . . . . . . . . . . . . .

5 5 9

4 Diferenciální metody 4.1 Bolzanova metoda . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Newtonova metoda . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Metoda regula falsi . . . . . . . . . . . . . . . . . . . . . . . .

14 14 17 20

5 Závěr

23

1

1

Úvod

Hledání minima funkce je často řešený problém optimalizace, oboru zabývajícího se určením nejlepšího řešení určitého matematicky popsatelného problému. I když v některých případech chceme znát maximum funkce, vztah max(f (x)) = −min(−f (x)) ukazuje, že jsou tyto pojmy spolu provázané. Nezáleží na tom, jestli potřebujeme znát maximum, nebo minimum, jedno můžeme určit z druhého. Optimalizační úlohy se vyskytují všude tam, kde je možné vybírat si z více možných rozhodnutí, a přitom kvalitu jednotlivých rozhodnutí ohodnotit nějakým reálným číslem. Konkrétní oblasti využití jsou: Matematika (teorie aproximace, optimalizace numerických procesů,. . .), Geometrie (geodézie, minimální křivky a plochy, optimální oblasti a tvary,. . .), Ekonomické a politologické teorie (využívání zdrojů a zásob, optimální skladba výroby, tvorba cen, rozložení rizika, strategické hry, teorie eskalace konfliktu,. . .), Fyzika (mechanika, geometrická optika, teorie pružnosti, hydrodynamika, teorie relativity,. . .), Přírodní vědy (modely fyzikálních, chemických a biologických procesů,. . .), Teorie řízení (optimální řízení, optimální systémy, hierarchické řízení, koordinační strategie,. . .), Teorie konstruování (optimalizace konstrukcí, optimalizace tvarů, optimální odhad neznámých parametrů, optimalizace dynamických vlastností mechanických systému, optimalizace spolehlivosti a rizika konstrukcí,. . .) a další jako například v dopravě, logistice a kdekoliv jinde, kde se nám podaří matematicky zformulovat optimalizační úlohu. Není vždy možné určit minimum přesně, nebo je to někdy příliš pracné, nám přitom stačí „přibližnáÿ hodnota, proto hledáme minimální hodnoty funkcí numericky. My se omezíme pouze na funkce jedné proměnné. V textu jsou uvedena tvrzení a věty, které v mé práci nejsou dokázány. Důvodem je skutečnost, že tato tvrzení jsou všeobecně známá (tudíž snáze uvěřitelná), případně by jejich dokazování „zbytečně natahovaloÿ text, v mnohých případech vedlo k definicím nových pojmů a uvádění tvrzení s nimi spjatých, což by změnilo podstatu textu, v některých případech zřejmě i odradilo nejednoho čtenáře. U každého tvrzení (a věty) jsou uváděny odkazy na příslušnou literaturu, v níž si zájemci mohou jednotlivé důkazy dohledat.

2

2

Pojmy a předpoklady

Definice 1. O funkci f řekneme, že má v bodě x = p lokální minimum, jestliže existuje otevřený interval I obsahující bod p takový, že f (p) ≤ f (x) pro všechna x ∈ I. Jestliže navíc f (p) < f (x) pro všechna x ∈ I, pak tomuto lokálnímu minimu říkáme ostré. Definice 2. Nechť je funkce f definována na intervalu I. (i) Jestliže pro každé dva body x, y ∈ I takové, že x < y, platí f (x) < f (y), nazývá se f rostoucí na intervalu I. (ii) Jestliže pro každé dva body x, y ∈ I takové, že x < y, platí f (x) > f (y), nazývá se f klesající na intervalu I. Definice 3. Funkce f se nazývá unimodální na intervalu I, jestliže má na tomto intervalu právě jeden lokální extrém. Zadanou funkci, jejíž minimum hledáme, budeme označovat jako účelovou. V celé práci předpokládáme, že účelová funkce je na zadaném intervalu spojitá a unimodální, přičemž nabývá (ostrého) lokálního minima v jeho vnitřním bodě, čili má na daném intervalu tvar přibližně odpovídající obrázku 2.1.

Obr. 2.1: Unimodální funkce

3

Jestliže funkce f nabývá své minimální hodnoty v bodě p a my tuto hodnotu hledáme na intervalu [a, b], pak podle výše uvedených předpokládů, je funkce f na intervalu [a, p] klesající a na intervalu [p, b] rostoucí. Věta 1. Nechť funkce f má v bodě p lokální minimum a nechť existuje derivace f ′ (p). Pak f ′ (p) = 0. [3, str. 187] Body, v nichž derivace funkce nabývá nulové hodnoty, označujeme jako stacionární. Stacionární bod nemusí být vždy lokálním extrémem (minimem nebo maximem), vzhledem k našim předpokladům však stacionární bod bude vždy dokonce lokálním minimem funkce f na intervalu [a, b]. Numerické metody hledání minima funkce jedné proměnné rozdělujeme na přímé (bez dalších požadavků) a metody diferenciální (derivační), v nichž převádíme problém hledání minima funkce na problém nalezení kořene derivace účelové funkce, čímž tedy navíc požadujeme existenci derivace na zadaném intervalu. Věta 2. Předpokládejme, že funkce f je spojitá na intervalu I = [a, b] a je na něm diferencovatelná. (i) Jestliže f ′ (x) > 0 pro všechna x ∈ I, pak f je rostoucí na I. (ii) Jestliže f ′ (x) < 0 pro všechna x ∈ I, pak f je klesající na I. [3, str. 183] V závěrečném srovnání metod pro nás bude užitečná definice řádu metody, tedy „míryÿ rychlosti její konvergence. Metody, které zde uvedeme, jsou vícekrokové. V každém kroku (iteraci) výpočtu je přitom spočítána aproximace hledaného řešení. Pro jednoduchost budou všechny funkce funkcemi proměnné x, aproximaci minima v k-tém kroku budeme značit xk , přesné řešení pak x∗ . Uveďme dále, že v některých literaturách (jako např. v [1]) je k-tá aproximace značena xk . Definice 4. Chybu ek k-té iterace definujeme vztahem ek = |xk − x∗ |. Předpokládejme nyní, že metoda je konvergentní, tzn. limk→∞ xk = x∗ . Existuje-li reálné číslo p ≥ 1 tak, že platí limk→∞

ek+1 epk

= C 6= 0,

řekneme, že daná iterační metoda je řádu p. Konstanta C se nazývá asymptotickou konstantou chyby.

4

3 3.1

Přímé metody Zlatý řez

Metoda zlatého řezu patří mezi postupné (adaptivní) komparativní metody hledání minima libovolné spojité unimodální funkce. Postupné komparativní metody spočívají v utvoření takové posloupnosti intervalů [ak , bk ], x∗ ∈ Ik = [ak , bk ], kde x∗ je hledané minimum, která vyhovuje vztahu: Ik ⊂ Ik−1 ⊂ Ik−2 ⊂ · · · ⊂ I2 ⊂ I1 ⊂ I = [a, b] Na počátku je zadán interval [a, b]. Každý následující interval je subintervalem předcházejícího, a proto posloupnost {Ik } konverguje. Metoda zlatého řezu spočívá v rozdělení každého intervalu tak, aby poměr větší části k menší byl roven poměru celého děleného intervalu k větší části. Vyznačme si zlatý řez na úsečce a = AB:

Obr. 3.1: Zlatý řez Za předpokladu, že délka úsečky a je rovna 1, odvodíme: r a−r

=

a r

r 1−r

=

1 r

r2 + r − 1 = 0 Kladný kořen značí délku úsečky r (přibližně 0,618).

5

Interval [a, b], v němž minimum hledáme, rozdělíme dvěma body c, d tak, aby platilo: c = a + (1 − r)(b − a)

a

d = a + r(b − a)

Jestliže f (c) ≤ f (d), pak se minimum nachází v subintervalu [a, d], s nímž provedeme následující iteraci výpočtu, tedy [a1 , b1 ] = [a, d]. V případě, že f (c) > f (d), minimum se nachází v subintervalu [c, b], tedy [a1 , b1 ] = [c, b], viz schematické obrázky 3.2 a 3.3. Vybarvená část plochy pod grafem funkce f značí výběr intervalu pro další výpočet.

Obr. 3.2: f (c) ≤ f (d)

Obr. 3.3: f (c) > f (d)

V průběhu výpočtu se zadaný interval I = [a, b] zmenšuje. Označme jeho velikost l. Pak platí: |[a, b]| = l = b − a |[a1 , b1 ]| = l1 = b1 − a1 = r · (b − a) |[a2 , b2 ]| = l2 = b2 − a2 = r 2 · (b − a) .. . |[ak , bk ]| = lk = bk − ak = r k · (b − a)

Pro odchylku (též chybu) výpočtu platí ε ≥ intervalu po poslední provedené iteraci.

6

1 2

· ln , kde ln je velikost

Z předchozích vztahů odvodíme počet kroků výpočtu: ≥ 12 · ln = ≥ rn .. .

ε 2ε b−a ) log( b−a 2ε log( 1r )

1 2

· r n · (b − a)

≤n

Pro dosažení přesnosti ε musíme provést aspoň N kroků výpočtu, kde N≥

log( b−a ) 2ε . − log r

Výsledek výpočtu (značíme x∗ ) zapíšeme s pomocí bodu xn , což je střed intervalu [an , bn ]: x∗ = xn ± 12 · ln

(případně x∗ = xn ± ε)

Příklad: Nalezněte minimum funkce f (x) = x + x32 s přesností ε = 0, 05, víte-li, že se nalézá v intervalu [0.5, 3].

Řešení: Nejdříve zjistíme počet kroků iterace nutných k zajištění přesnosti ε:

N≥

) log( b−a 2ε − log r

=

log( 2,5 ) 0,1 − log r

. = 6, 69 ⇒ N = 7

Když víme, že 7 iterací stačí k dosažení požadované přesnosti, nezbývá, než je provést a určit přibližné řešení. 1.krok: I = [0.500000, 3.000000], c = 1, 454915, d = 2, 045085, f (c) = 2, 872163, f (d) = 2, 762381 f (c) > f (d) ⇒ I1 = [1.454915, 3.000000] 7

2.krok: I1 = [1.454915, 3.0000000], c = 2, 045085, d = 2, 409830, f (c) = 2, 762381, f (d) = 2, 962423 f (c) ≤ f (d) ⇒ I2 = [1.454915, 2.409830] 3.krok: I2 = [1.454915, 2.409830], c = 1, 819660, d = 2, 045085, f (c) = 2, 725686, f (d) = 2, 762381 f (c) ≤ f (d) ⇒ I3 = [1.454915, 2.045085] 4.krok: I3 = [1.454915, 2.045085], c = 1, 680340, d = 1, 819660, f (c) = 2, 742835, f (d) = 2, 725686 f (c) > f (d) ⇒ I4 = [1.680340, 2.045085] 5.krok: I4 = [1.680340, 2.045085], c = 1, 819660, d = 1, 905765, f (c) = 2, 725686, f (d) = 2, 731770 f (c) ≤ f (d) ⇒ I5 = [1.680340, 1.905765] 6.krok: I5 = [1.680340, 1.905765], c = 1, 766445, d = 1, 819660, f (c) = 2, 727882, f (d) = 2, 725686 f (c) > f (d) ⇒ I6 = [1.766445, 1.905765] 7.krok: I6 = [1.766445, 1.905765], c = 1, 819660, d = 1, 825549, f (c) = 2, 725686, f (d) = 2, 726691 f (c) ≤ f (d) ⇒ I7 = [1.766445, 1.852549] I7 = [1.766445, 1.852549] ⇒ x7 = 1, 809497 Přibližné řešení: x∗ = 1, 809497 ± 0, 043052 Zápis výsledku s požadovanou přesností: x∗ = 1, 81 ± 0, 05 Přesné řešení: x∗ =

√ . 3 6 = 1, 817120

Na obrázku 3.4 je graficky znázorněn předchozí výpočet. Zmenšující se křivky pod grafem účelové funkce (jež kopírují) udávají prohledávaný interval v jednotlivých iteracích.

8

Obr. 3.4: Metoda zlatého řezu pro f (x) = x +

3.2

3 . x2

Kvadratická interpolace

Metoda kvadratické interpolace spočívá v aproximaci účelové funkce na daném intervalu parabolou. K proložení funkce parabolou na daném intervalu potřebujeme znát funkční hodnoty ve třech bodech, které ji určí jednoznačně. Označme je r, s, t. Body r a t budou okrajové body daného intervalu, bod s bude ležet mezi nimi, přitom musí platit: f (r) > f (s) < f (t). Parabolu, jíž aproximujeme zadanou funkci, vyjádřeme Lagrangeovým polynomem. Q(x) =

f (r)(x−s)(x−t) (r−s)(r−t)

+

f (s)(x−r)(x−t) (s−r)(s−t)

+

f (t)(x−r)(x−s) (t−r)(t−s)

Parabola Q(x) má minimum v bodě x, pro nějž Q′ (x) = 0. Q′ (x) =

f (r)(2x−s−t) (r−s)(r−t)

+

f (s)(2x−r−t) (s−r)(s−t)

9

+

f (t)(2x−r−s) (t−r)(t−s)

Označme minimum paraboly m. Pak platí: m=

1 2

·

f (r)(t2 −s2 )+f (s)(r 2 −t2 )+f (t)(s2 −r 2 ) f (r)(t−s)+f (s)(r−t)+f (t)(s−r)

Grafické provedení předchozího postupu:

Obr. 3.5: Kvadratická interpolace

Předcházející postup může mít 3 různé výsledky (závisející na volbě bodů r, s, t) lišících se polohou bodu m: 1) m < s: Pokud f (m) < f (s), pak se minimum nachází v intervalu [r, s], tudíž body r, s budou nové okrajové body intervalu r, t, bod s bude spočítaný bod m. Pokud f (m) > f (s), pak se minimum nachází v intervalu [m, t], tudíž body m, t budou nové okrajové body intervalu r, t, bod s „zůstane na místěÿ. 2) m = s: Algoritmus nás „dovedl ke konciÿ (tzn. lepší aproximaci tímto postupem se stejnými počátečními hodnotami nezískáme), bod s je minimem paraboly. Buď jsme jej zvolili tak nešťastně, že i když je minimem kvadratické funkce proložené body r,s,t, není minimem účelové funkce (zjistíme například zkoumáním funkčních hodnot v ε-okolí bodu s). Pak nám nezbývá, než zvolit jiný bod s a algoritmus opakovat. V opačném případě je minimum účelové funkce shodné s minimem paraboly, tudíž výpočet ukončíme, jelikož jsme našli přesnou hodnotu minima.

10

3) m > s: Pokud f (m) < f (s), pak se minimum nachází v intervalu [s, t], tudíž body s, t budou nové okrajové body intervalu r, t, bod s bude spočítaný bod m. Pokud f (m) > f (s), pak se minimum nachází v intervalu [r, m], tudíž body r, m budou nové okrajové body intervalu r, t, bod s „zůstane na místěÿ. U metody kvadratické interpolace neurčujeme počet kroků nutných k dosažení určité chyby (přesnosti). Výpočet končí ve chvíli, kdy |xk −xk−1 | ≤ ε, kde xn je minimum m paraboly Q(x) v n-tém kroku výpočtu. Výsledek zapíšeme ve tvaru x∗ ≈ xk .


Řešení: 1.krok: [a, b] = [0.500000, 3.000000], r = 0.500000, s = 1.500000, t = 3.000000, x0 = 1.750000 f (r) = 12.500000, f (s) = 2.833333, f (t) = 3.333333 m = 2.208333, f (m) = 2.823499 (s < m) ∧ (f (m) < f (s)) ⇒ r = 1.500000, s = 2.208333, t = 3.000000 x1 = 2.208333, |x1 − x0 | = 0.458333 2.krok: [a, b] = [1.500000, 3.000000], r = 1.500000, s = 2.208333, t = 3.000000, x1 = 2.208333 f (r) = 2.833333, f (s) = 2.823499, f (t) = 3.333333 m = 1.869995, f (m) = 2.727902 (m < s) ∧ (f (m) < f (s)) ⇒ r = 1.500000, s = 1.869995, t = 2.208333 x2 = 1.869995, |x2 − x1 | = 0.338339

11

3.krok: [a, b] = [1.500000, 2.208333], r = 1.500000, s = 1.869995, t = 2.208333, x2 = 1.869995 f (r) = 2.833333, f (s) = 2.727902, f (t) = 2.823499 m = 1.862831, f (m) = 2.727350 (m < s) ∧ (f (m) < f (s)) ⇒ r = 1.500000, s = 1.862831, t = 1.869995 x3 = 1.862831, |x3 − x2 | = 0.007163 Přibližné řešení: x∗ ≈ 1, 862831 Zápis výsledku s požadovanou přesností: x∗ = 1, 86 ± 0, 05 Přesné řešení: x∗ =

√ . 3 6 = 1, 817120

Na obrázku 3.6 je graficky znázorněn předchozí výpočet. Jednotlivé aproximace jsou znázorněny postupně se zvětšujícími „čárkamiÿ protínajícími graf účelové funkce, minimum je znázorněno nejdelší „čárkouÿ.

Obr. 3.6: Metoda kvadratické interpolace pro f (x) = x +

12

3 . x2

Jak vidíme z výpočtu i z obrázku, metoda nemusí dávat přesné (tzn. pravdivé) řešení. V našem případě je přesnost ε zvolena vhodně, tudíž zápis výsledku s požadovanou přesností je pravdivý. Uvážíme-li však například přesnost ε = 0, 01, výpočet bude stejný, ovšem zápis výsledku s požadovanou přesností nebude pravdivý. Na vině je podmínka ukončení výpočtu, v níž požadujeme, aby rozdíl dvou následujících aproximací byl menší než daná chyba. To ovšem nezaručuje, že jsme též „dostatečně blízkoÿ skutečného minima účelové funkce. K zaručení skutečné přesnosti (pravdivosti) výsledku je třeba brát přesnost jako interval [r, t], který se postupně zmenšuje a kde máme jistotu, že obsahuje skutečné minimum účelové funkce. To však vede k podstatně delšímu výpočtu (v našem případě by to znamenalo 27 iterací). Později uvidíme, že můžeme k odhadu chyby navíc použít i jiné vztahy, vzhledem k našim předpokladům (kdy nepožadujeme existenci derix −x vace účelové funkce) z nich můžeme využít pouze | k+1xk k | < ε, který stejně pravdivost výsledku nemusí zaručit (viz. náš příklad a ε = 0, 01). Závěrem dodejme, že existuje více implementací metody kvadratické interpolace, které se liší především volbou vnitřního bodu s intervalu [r, t]. Známá je pak především implementace s tzv. ekvidistantními body, kde bod s je vždy středem intervalu [r, t]. To dovoluje jednodušší výpočet hodnoty minima m paraboly Q(x), k dosažení stejné přesnosti je však potřeba více kroků výpočtu. Navíc se musí zajistit, aby nenastala situace f (r) = f (t), pak by totiž v následující iteraci bod s splynul s bodem m.

13

4

Diferenciální metody

Druhou skupinou metod k hledání minima funkce jedné proměnné jsou metody derivační, nebo-li diferenciální. Doposud jsme po zadané účelové funkci požadovali „pouzeÿ, aby na daném intervalu byla unimodální a měla tvar odpovídající obrázku 2.1. V derivačních metodách budeme navíc požadovat, aby funkce byla na zadaném intervalu diferencovatelná, tedy aby byla její derivace definovaná v každém bodě tohoto intervalu.

4.1

Bolzanova metoda

Nejjednodušší iterační metodou pro hledání minima je Bolzanova metoda, též známá jako metoda bisekce čili půlení intervalů. Jak název napovídá, prohledávaný interval se s každým krokem výpočtu zmenší na polovinu v závislosti na hodnotě derivace funkce, jejíž minimum hledáme. Na schematickém obrázku 4.1 je znázorněna obecná funkce f na intervalu [a, b] spolu se svojí derivací (dole) na tomtéž intervalu. Dále je na obrázku naznačeno půlení zadaného intervalu v prvních třech krocích. Jestliže označíme p bod, v němž funkce f nabývá svého minima na intervalu [a, b], pak je f (jak již bylo zmíněno v oddílu 2) klesající na intervalu [a, p] a rostoucí na intervalu [p, b]. Podle věty 2 je tedy ∀x ∈ [a, p] : f ′ (x) < 0, ∀x ∈ [p, b] : f ′ (x) > 0 a podle věty 1 v bodě p je f ′ (x) = 0. Položme a0 = a, b0 = b, x0 = 12 (a0 + b0 ). Pro k = 0, 1, 2, ... opakujeme následující postup, dokud nedocílíme požadované přesnosti: (i) f ′ (xk ) < 0 ⇒ p > xk ⇒ ak+1 = xk , bk+1 = bk , xk+1 = 21 (ak+1 + bk+1 ) (ii) f ′ (xk ) > 0 ⇒ p < xk ⇒ ak+1 = ak , bk+1 = xk , xk+1 = 21 (ak+1 + bk+1 ) (iii)f ′ (xk ) = 0 ⇒ xk je minimum Jestliže označíme lk = |[ak , bk ]|, pak pro chybu výpočtu platí: ε ≥ Po N-tém kroku výpočtu: x∗ = xN ± 21 · lN , kde lN = 2−N (b − a). Nutný počet iterací k dosažení přesnosti ε: N ≥ 14

log b−a 2ε log 2

.

1 2

· lk .

Obr. 4.1: Bolzanova metoda


Řešení: Nejdříve zjistíme počet kroků iterace nutných k zajištění přesnosti ε: N≥

log b−a 2ε log 2

f ′ (x) = 1 −

=

2,5 log( 0,1 )

log 2

. = 4, 64 ⇒ N = 5

6 x3

Víme, že 5 kroků výpočtu stačí k zajištění přesnosti ε, víme, jak určit funkční hodnoty derivace funkce f , nic nám nebrání aplikovat výše uvedený postup. 1.krok: I = [0.500000, 3.000000], a0 = 0.500000, b0 = 3.000000, x0 = 1.750000, f ′ (x0 ) = −0.119534 < 0 ⇒ I1 = [1.750000, 3.000000] 15

2.krok: I1 = [1.750000, 3.000000], a1 = 1.750000, b1 = 3.000000, x1 = 2.375000, f ′(x1 ) = 0.552121 > 0 ⇒ I2 = [1.750000, 2.375000] 3.krok: I2 = [1.750000, 2.375000], a2 = 1.750000, b2 = 2.375000, x2 = 2.062500, f ′(x2 ) = 0.316137 > 0 ⇒ I3 = [1.750000, 2.062500] 4.krok: I3 = [1.750000, 2.062500], a3 = 1.750000, b3 = 2.062500, x3 = 1.906250, f ′(x3 ) = 0.133813 > 0 ⇒ I4 = [1.750000, 1.906250] 5.krok: I4 = [1.750000, 1.906250], a4 = 1.750000, b4 = 1.906250, x4 = 1.828125, f ′(x4 ) = 0.017950 > 0 ⇒ I5 = [1.750000, 1.828125] Přibližné řešení: x∗ = x5 ± 0, 039063 = 1, 789063 ± 0, 039063 Zápis výsledku s požadovanou přesností: x∗ = 1, 79 ± 0, 05 Přesné řešení: x∗ =

√ . 3 6 = 1, 817120

Na obrázku 4.2 je graficky znázorněn předchozí výpočet. Zmenšující se křivky pod grafem účelové funkce (jež kopírují) udávají prohledávaný interval v jednotlivých iteracích.

Obr. 4.2: Bolzanova metoda pro f (x) = x +

16

3 . x2

4.2

Newtonova metoda

Další iterační metodou, která hledá kořen funkce (v našem případě kořen derivace), je Newtonova metoda (též Newton-Raphsonova metoda, či metoda tečen). Metoda využívá navíc i druhé derivace, tudíž musíme zvýšit naše požadavky na účelovou funkci o existenci její druhé derivace v každém bodě zadaného intervalu [a, b]. Jeden z názvů napovídá, že metoda využívá k nalezení kořenu tečen, jak je mimo jiné vidět i z geometrické interpretace naznačené na obrázku 4.3. Vyznačena je opět obecná funkce f a pod ní její derivace na zadaném intervalu [a, b], navíc jsou na obrázku vidět tři po sobě jdoucí iterace v grafickém provedení.

Obr. 4.3: Newtonova metoda

17

Křivku v okolí kořene nahrazujeme tečnou v bodě [xi , f (xi )]. Průsečík tečny s osou x je nová aproximace xi+1 . Vycházíme z toho, že z hodnoty první derivace funkce v určitém bodě se dá určit rovnice tečny procházející tímto bodem ([3, str. 159]). V našem případě pak odvodíme následující vztah: −f ′ (xk ) xk+1 −xk

= f ′′ (xk ) ⇒ xk+1 = xk −

f ′ (xk ) f ′′ (xk )

U metody tečen neurčujeme počet kroků výpočtu. Výpočet končí ve chvíli, kdy |xk+1 − xk | ≤ ε, kde ε je požadovaná přesnost. Značnou nevýhodou metody je fakt, že nekonverguje vždy. Postačující podmínky konvergence ([5], [1, str. 44]): (i) f ′ (a)f ′ (b) < 0 (ii) ∀x ∈ [a, b] : sgn(f ′′(x)) = konst (6= 0) (iii)∀x ∈ [a, b] : sgn(f ′′′(x)) = konst (6= 0) Za počáteční aproximaci je vhodné volit bod x0 tak, že : f ′ (x0 )f ′′′ (x0 ) > 0. Při nevhodné volbě počátečního bodu může další aproximace ležet mimo interval [a, b], či blíže k předcházející aproximaci než je požadovaná přesnost, avšak „dostatečněÿ daleko od kořene (na to, abychom dostali špatný výsledek). Konvergence Newtonovy metody je tím rychlejší, čím více se účelová funkce blíží kvadratické parabole, pro niž dostaneme přesné řešení po jednom kroku. Příklad: Nalezněte minimum funkce f (x) = x + x32 s přesností ε = 0, 05, víte-li, že se nalézá v intervalu [0.5, 3].

Řešení: f ′ (x) = 1 −

6 , x3

f ′′ (x) =

18 x4

Zvolme x0 = 1.750000 (1. aproximace z metody bisekce) 1.krok: x0 = 1, 750000, f ′(x0 ) = −0, 119534, f ′′(x0 ) = −1, 919200 ⇒ ⇒ x1 = 1, 812283, |x1 − x0 | = 0, 062283 2.krok: x1 = 1, 812283, f ′(x1 ) = −0, 008029, f ′′(x1 ) = 1, 668662 ⇒ ⇒ x2 = 1, 817095, |x2 − x1 | = 0, 004812 18

Přibližné řešení: x∗ ≈ 1, 817095 Zápis výsledku s požadovanou přesností: x∗ = 1, 82 ± 0, 05 Přesné řešení: x∗ =

√ 3

. 6 = 1, 817120


Obr. 4.4: Newtonova metoda pro f (x) = x +

3 . x2

Podobně jako u metody kvadratické interpolace nemáme ani zde zaručenu požadovanou přesnost výpočtu, jelikož výpočet opět končí ve chvíli, kdy jsou si „dostatečněÿ blízké dvě následující aproximace. Možných chyb se můžeme vyvarovat například stanovením větší přesnosti, či „dostatečněÿ blízkou aproximací skutečnému minimu. Zisk počáteční aproximace (počátečních aproximací) se zpravidla provádí metodou bisekce, vzhledem k rychlosti konvergence (více v sekci Závěr) Newtonovy metody nemá zvětšení přesnosti velký vliv na délku výpočtu. Dále je možné pro odhad chyby použít i vztahy | xk+1xk−xk | < ε, |f ′(xk )| < ε. 19

4.3

Metoda regula falsi

V předchozí metodě potřebujeme v každém kroku spočítat jak první, tak druhou derivaci funkce v daném bodě. Jelikož výpočet derivace funkce nemusí být vždy snadný, nahradíme derivaci aproximací: f ′ (xk ) ≈

f (xk )−f (xk−1 ) xk −xk−1

V našem případě však potřebujeme znát kromě první derivace i druhou, tudíž stejnou aproximaci provedeme „ažÿ pro druhou derivaci, první bude zastávat roli „původníÿ funkce ve výše uvedené aproximaci. Dostaneme tedy: f ′′ (xk ) ≈

f ′ (xk )−f ′ (xk−1 ) xk −xk−1

Nahradíme-li v iteračním vzorci Newtonovy metody druhou derivaci uvedenou aproximací, dostáváme: xk+1 = xk −

xk −xk−1 f ′ (xk ) f ′ (xk )−f ′ (xk−1 )

Právě uvedený vzorec reprezentuje tzv. metodu sečen. Z názvu plyne její geometrický význam. Podobně jako Newtonova metoda není metoda sečen vždy konvergentní ([1, str. 49]). Chceme-li zajistit konvergenci metody pro libovolné počáteční aproximace (všimněme si, že narozdíl od Newtonovy metody nyní potřebujeme znát dvě předchozí aproximace k výpočtu té následující), je třeba, aby platilo: f ′ (xk )f ′ (xk−1 ) < 0. Poté metoda s předpisem xk+1 = xk −

xk −xs f ′ (xk ), f ′ (xk )−f ′ (xs )

kde s = s(k) je největší index takový, že f ′ (xk )f ′ (xs ) < 0, je konvergentní ∀x ∈ [a, b], přičemž předpokládáme, že počáteční aproximace x0 , x1 jsou vybrány tak, že f ′ (x0 )f ′ (x1 ) < 0 (tj. např. x0 = a, x1 = b). [1, str. 51] Uvedená metoda se nazývá metoda regula falsi. Její geometrický význam je velmi podobný metodě sečen. Aproximaci xk+1 získáme jako průsečík přímky vedené body [xk , f ′ (xk )], [xs , f ′ (xs )] s osou x. Na rozdíl od metody sečen aplikujeme další postup na ten z intervalů [xk , xk+1 ], [xs , xk+1 ], v jehož 20

koncových bodech má funkce f ′ opačná znaménka (pokud f ′ (xk+1 ) = 0, xk+1 je hledané minimum funkce f ). Na obrázku 4.5 je schematicky zobrazen popsaný postup.

Obr. 4.5: Metoda regula falsi Podobně jako u Newtonovy metody výpočet ukončíme ve chvíli, kdy |xk+1 − xk | ≤ ε, kde ε je požadovaná přesnost. Dodejme, že pro odhad chyby x −x se používají i vztahy | k+1xk k | < ε, |f ′ (xk )| < ε. Příklad: Nalezněte minimum funkce f (x) = x + x32 s přesností ε = 0, 05, víte-li, že se nalézá v intervalu [0.5, 3].

Řešení: Pro nalezení prvních dvou iterací můžeme použít např. metodu bisekce. Tedy: x0 = 1, 750000, x1 = 2, 375000. Ověříme,že f ′ (x0 )f ′ (x1 ) < 0. f ′ (x0 ) = f ′ (1, 750000) = −0, 119534, f ′(x1 ) = f ′ (2, 375000) = 0, 552121 ⇒ ⇒ f ′ (x0 )f ′ (x1 ) < 0. 1.krok: x2 = x1 −

x1 −x0 f ′ (x1 ) f ′ (x1 )−f ′ (x0 )

= 1, 861230, |x2 − x1 | = 0, 062283

f ′ (x0 ) = f ′ (1, 750000) = −0, 119534, f ′(x2 ) = f ′ (1, 861230) = 0, 069426 ⇒ ⇒ f ′ (x0 )f ′ (x2 ) < 0. 2.krok: x3 = x1 −

x2 −x0 f ′ (x2 ) f ′ (x2 )−f ′ (x0 )

= 1, 820363, |x3 − x2 | = 0, 040867, f ′ (x3 ) = 0, 005334

21

Přibližné řešení: x∗ ≈ 1, 820363 Zápis výsledku s požadovanou přesností: x∗ = 1, 82 ± 0, 05 Přesné řešení: x∗ =

√ 3

. 6 = 1, 817120


Obr. 4.6: Metoda regula falsi pro f (x) = x +

3 . x2

V našem příkladu jsme metodu ukončili ve chvíli, kdy |xk+1 − xk | ≤ ε, navíc jsme otestovali, že |f ′(xk+1 )| < ε. Jelikož může dojít i u této metody k tomu, že jsou si aproximace xk , xk+1 bližší než ε, přičemž vzdálenost od skutečného minima je větší, pravidelně se užívá při testování konce výpočtu současně více vztahů uvedených před příkladem.

22

5

Závěr

V práci bylo popsáno pět numerických metod pro nalezení minima funkce jedné proměnné. Při jejich porovnávání nás přirozeně budou zajímat kritéria: i) Jaké jsou požadavky na tu kterou metodu? ii) Pro která x ze zadaného intervalu [a, b] metoda konverguje ke skutečnému minimu funkce f a jak rychle? Odpovědi na první otázku už známe. V první kapitole jsme si řekli, že budeme požadovat, aby účelová funkce měla daný tvar (jediný ostrý lokální extrém, a to minimum nacházející se v intervalu (a, b)). Tento požadavek byl stejný pro všechny metody. Základním rozdílem, kterým jsme metody rozdělili do dvou skupin, byl požadavek na derivaci účelové funkce. Jak víme, přímé metody ke svému „choduÿ nepotřebovali žádnou znalost o derivaci účelové funcke, zatímco metody diferenciální ze znalosti (pochopitelně také z existence) derivace účelové funkce vycházely, u Newtonovy metody se přitom jednalo dokonce i o druhou derivaci účelové funkce. Z popsaného tedy vyplývá zřejmá výhoda prvně zmíněných (přímých) metod. Druhé kritérium ještě zcela zodpovězeno nemáme. Víme již, že metody zlatého řezu, bisekce (Bolzanova metoda) a regula falsi konvergují vždy k minimu účelové funkce, přičemž u metody regula falsi musí první dvě aproximace x0 , x1 splňovat podmínku f ′ (x0 )f ′ (x1 ) < 0, což můžeme vždy zajistit. Na Newtonovu metodu a metodu kvadratické interpolace již musíme klást další nároky, aby konvergovaly k minimu účelové funkce. U prvně zmíněné metody se jedná o zvolení „dostatečněÿ blízké aproximace ke skutečnému minimu, což zpravidla provádíme metodou bisekce. Konvergence metody kvadratické interpolace pak závisí na volbě bodů r, s, t, nejvíce pak na volbě bodu s. Nebýt rychlosti konvergence jistě bychom za nejužitečnější zvolili metodu zlatého řezu, jež vždy konverguje, a přitom má nejméně požadavků, zatímco za nejméně užitečnou bychom patrně považovali Newtonovu metodu, která má nejvíce požadavků, přičemž ani nekonverguje vždy. Právě rychlost konvergence však ukazuje podstatnou výhodu Newtonovy metody, která je jako jediná z uvedených metodou druhého řádu (viz druhá kapitola, též je možné setkat se s pojmem kvadratická konvergence) [1, str. 41]. Metody zlatého řezu, bisekce a regula falsi jsou metodami prvního řádu (též lineárně konvergentní) [1, str. 51], problém konvergence metody kvadratické interpolace je dosti složitý. Z uvedeného tedy plyne, že ač má Newtonova metoda nejvíce 23

požadavků, konverguje „značněÿ rychleji než ostatní metody. Vybrat některou z metod jako nejvhodnější pro všechny funkce jedné proměnné je nemožné. Bylo by reálné se pro jednu pevně rozhodnout na základě určitých potřeb (vždy konvergentní, rychlost výpočtu,. . .), celkově má však každá metoda nějaké klady a nějaké zápory. Proto bych doporučil zvolit metodu až podle předpisu funkce, jejíž minimum chceme najít, případně (a to je zřejmě nejvýhodnější) užít k nalezení minima kombinaci dvou metod jako například Bolzanovy a Newtonovy metody, či metod zlatého řezu a kvadratické interpolace apod. Na přiloženém CD jsou kromě souborů s textem práce i programy psané v Matlabu demonstrující jednotlivé metody na několika funkcích.

24

Literatura

[1] Horová Ivana, Zelinka Jiří. Numerické metody. 2.vyd. Masarykova univerzita Brno, 2004. 294s. ISBN 80-210-3317-7. [2] Mathews, John H. Numerical methods for mathematics, science and engineering. 2.vyd. Englewood Cliffs: Prentice-Hall International, 1992. 646s. X.ISBN.0-13-625047-5. [3] Novák, Vítězslav. Diferenciální počet v R. 2.vyd. Masarykova univerzita Brno, 1997. 231s. ISBN 80-210-1561-6. [4] Svobodová Vařeková Radka. Optimalizace. Masarykova univerzita Brno. http://ncbr.chemi.muni.cz/˜ n19n/vyuka/optimalizace/ [5] Vítečková Miluše, Jedlička David. Optimalizace funkce jedné proměnné. Vysoká škola báňská - Technická univerzita Ostrava, 2003. http://www.fs.vsb.cz/books/StatickaOptimalizace/2fjedne/2fjedne2.htm

25

Hledání minima funkce jedné proměnné

Recommend Documents