UNIVERZITA PALACKÉHO V OLOMOUCI PŘÍRODOVĚDECKÁ FAKULTA KATEDRA MATEMATICKÉ ANALÝZY A APLIKACÍ MATEMATIKY
BAKALÁŘSKÁ PRÁCE Numerické metody jednorozměrné minimalizace
Vedoucí bakalářské práce: RNDr. Horymír Netuka, Ph.D. Rok odevzdání: 2013
Vypracovala: Barbora Dohnalová MATAP, III. ročník
Prohlášení Prohlašuji, že jsem tuto práci zpracovala samostatně. Veškeré zdroje a literatura, kterých bylo při práci používáno, jsou citovány a uvedeny v seznamu literatury.
V Olomouci dne 26.4.2013
Poděkování Ráda bych na tomto místě poděkovala svému vedoucímu práce, RNDr. Horymíru Netukovi, PhD., za cenné připomínky a rady, bez nichž by práce nemohla být vypracována. Také děkuji své rodině a příteli, kteří mě během psaní podporovali.
Obsah Úvod
4
1 Základní pojmy
5
2 Typy metod
8
3 Metody druhého řádu 3.1 Newtonova metoda . . . . . . . . . . . . . . . . . . . . . . . . . .
10 10
4 Metody prvního řádu 4.1 Metoda půlení intervalu . . 4.2 Metoda sečen . . . . . . . . 4.3 Metoda regula falsi . . . . . 4.4 Metoda kubické interpolace
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
12 12 13 13 15
5 Metody nevyužívající derivace 5.1 Metoda kvadratické interpolace . . 5.2 Rovnoměrná komparativní metoda 5.3 Metoda zlatého řezu . . . . . . . . 5.4 Fibonacciho metoda . . . . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
18 18 21 23 25
. . . .
. . . .
. . . .
6 Příklady
26
Literatura
36
Úvod Tématem této práce jsou numerické metody, které se používají při hledání nepodmíněného minima funkcí jedné proměnné. Numerické metody řeší mnoho matematických úloh, u kterých neznáme přesné analytické řešení nebo je nalezení takového řešení těžce proveditelné. Hledání minima funkce je součástí procesu optimalizace, jenž má široké využití. Optimalizaci nalezneme například v přírodních vědách (optimalizace chemických či biologických procesů) nebo v ekonomice (optimalizace výrobních procesů, využití skladových zásob). Cílem práce je seznámení se nejenom se základními numerickými metodami, jenž jsou obsahem studia na Přírodovědecké fakultě Univerzity Palackého, jako je například Newtonova metoda či metoda sečen, ale také s méně známými numerickými metodami, jimž nejsou věnovány učební texty, jako je kvadratická či kubická interpolace. Záměrem této práce není srovnání daných metod ani zkoumání jejich rychlosti konvergence.
4
1. Základní pojmy V této kapitole si zavedeme základní pojmy vztahující se k tématu jednorozměrné minimalizace. Zároveň předpokládáme, že čtenář má základní znalosti z oblasti numerických metod a matematické analýzy, tedy chápe pojmy jako rostoucí funkce, derivace funkce, iterace apod. Definice 1. Řekneme, že funkce f má v bodě x0 nepodmíněné lokální minimum (resp. maximum), jestliže existuje δ-okolí bodu x0 , takové, že f (x0 ) ≤ f (x) (resp.f (x) ≥ f (x0 )) pro každé x z δ-okolí. Bude-li platit ostrá nerovnost, nazveme lokální minimum (resp. maximum) ostré. Definice 2. Nechť f je funkce s definičním oborem D(f ), M ⊂ D(f ) je neprázdná množina. Řekneme, že funkce f má v bodě x0 nepodmíněné globální minimum (resp. maximum) na M , jestliže ∀ ∈ M platí f (x) ≥ f (x0 ) (resp. f (x) ≤ f (x0 )). Bude-li platit ostrá nerovnost, nazveme globální minimum (resp. maximum) ostré. Definice 3. Funkce f se nazývá unimodální na intervalu I, jestliže má na I právě jeden lokální extrém. (viz [10], str.3) Poznámka 1. Unimodální funkce se dá definovat více způsoby, další z nich je ukázán v následující definici. Definice 4. Funkce f : R −→ R se nazývá unimodální na intervalu I, jestliže existuje bod x∗ ∈ I takový, že x∗ je bodem minima funkce f na I a jestliže pro libovolné x1 , x2 ∈ I takové, že x1 < x2 platí x2 ≤ x∗ =⇒ f (x1 ) > f (x2 ) x1 ≥ x∗ =⇒ f (x2 ) > f (x1 ). (viz [5], str.22)
5
Obrázek 1: Příklad unimodální funkce
Definice 5. Nechť je funkce f definovaná na intervalu I ⊆ R. Pokud pro libovolnou dvojici bodů x1 , x2 ∈ I a libovolné λ ∈ [0, 1] platí f (λx1 + (1 − λ)(x2 )) ≤ λf (x1 ) + (1 − λ)f (x2 ), pak funkci f nazýváme konvexní funkce na intervalu I. Poznámka 2. Dá se ukázat, že funkce, která je v intervalu [a, b] konvexní, je v tomto intervalu také unimodální. Účelovou funkcí nazýváme funkci, jejíž minimum hledáme.
Budeme se zabývat numerickými metodami pro unimodální spojité funkce, tedy budeme hledat řešení optimalizační úlohy x∗ = minf (x), x ∈ [a, b], kde f ∗ = f (x∗ ).
V teorii numerických metod je také třeba definovat přesnost vyznačení optimálního řešení x∗ , to jest číslo ε ≥ 0 takové, pro něž k-tá aproximace řešení xk vyhovuje nerovnosti |xk − x∗ | ≤ ε. Tento vztah bývá také nazýván odhadem absolutní chyby k-té aproximace. Přesností vyznačení optimální hodnoty účelové funkce f ∗ je číslo δ ≥ 0, pro něž 6
hodnota účelové funkce v bodě xk vyhovuje této nerovnosti |f (xk ) − f ∗ | ≤ δ.
Obrázek 2: Přesnost řešení
Splnění podmínek přesnosti je však značně problematické, jelikož je k jeho výpočtu třeba znát optimální řešení x∗ . Naším cílem je ale právě toto řešení najít, je tedy zřejmé, že při výpočtu numerických úloh pomocí algoritmů nelze dané vztahy pro přesnost užít. Z tohoto důvodu užíváme modifikované vztahy, které využívají hodnot xk , xk+1 , tedy hodnot, jichž dosáhneme výpočtem. Vztah užívající odhad absolutní chyby |xk+1 − xk | ≤ ε1 , a vztah pro odhad relativní chyby |xk+1 − xk | ≤ ε2 , |xk | kde ε1 ,ε2 jsou kladná reálná čísla.
7
2. Typy metod Numerické metody jednorozměrné minimalizace lze rozdělit na dva typy: • metody využívající derivace - metody vyžadující hodnoty účelové funkce a hodnoty první a druhé derivace této funkce ve všech bodech z intervalu [a, b] – metody druhého řádu - metody vyžadujicího hodnoty účelové funkce a hodnoty první i druhé derivace ∗ Newtonova metoda – metody prvního řádu - metody vyžadujicího hodnoty účelové funkce a hodnoty první derivace ∗ metoda půlení intervalu ∗ metoda sečen ∗ metoda regula falsi ∗ metoda kubické interpolace Metody využívající derivace jsou založeny na principu nalezení jediného stacionárního bodu, v němž má funkce f (x) na intervalu I ostré globální minimum, a který je reálným kořenem rovnice f 0 (x) = 0. V některých případech je rovnice řešitelná analyticky, jesliže však není algebraická, je třeba její řešení najít numerickými metodami. Pro zajištění nalezení minima je třeba předpokládat, že funkce f (x) je na intervalu I = [a, b] spojitá a f (a),f (b) mají rozdílná znaménka, tedy platí f (a)f (b) < 0.
8
• metody nevyužívající derivace - metody vyžadující jen hodnoty účelové funkce – metoda kvadratické interpolace – rovnoměrná komparativní metoda – metoda zlatého řezu – Fibonacciho metoda Metody nevyužívající derivace, někdy také nazývané jako metody nultého řádu, neužívají derivace minimalizovaných funkcí ani jejich aproximace. Byly využívany především v 60. letech 20. století, kdy byly považovány za spolehlivější než metody prvního a druhého řádu. Důvodem, proč je používáme dodnes, je fakt, že metody prvního a druhého řádu nefungují pro všechny optimalizační problémy a také to, že technická realizace přímých metod je snadná.
9
3. Metody druhého řádu 3.1. Newtonova metoda Newtonova metoda (neboli metoda tečen) je iterační metoda, ke které je třeba znát počáteční hodnotu xk a k ní jsme schopni vypočítat hodnoty f (xk ), f 0 (xk ), f 00 (xk ). Geometrický význam této metody spočívá v tom, že v počátečním bodě [xk , f (xk )] sestrojíme tečnu ke grafu funkce f (x). V bodě, ve kterém se protne graf tečny s osou x, leží bod [xk+1 , 0], protože platí (xk+1 − xk )f 00 (xk ) = −f 0 (xk ). Hledáme tedy průsečík přímky y s osou x 00
0
y = (xk+1 − xk )f (xk ) + f (xk ) = 0 => xk+1
f 0 (xk ) = xk − 00 . f (xk )
Z předpisu iteračního procesu je zřejmé, že nezávisí na hodnotě f (xk ). Newtonova metoda se tedy dá technicky zjednodušit na iterační řešení rovnic typu: g(x) = 0, kde g(x) ≡ f 0 (x). Poté se dá metoda psát stylem: xk+1 = xk −
g(xk ) . g 0 (xk )
Výpočet se ukončí při splnění nerovnosti |xk+1 − xk | ≤ ε.
10
Obrázek 3: Newtonova metoda
11
4. Metody prvního řádu 4.1. Metoda půlení intervalu Metoda půlení intervalu je známá také pod pojmy Bolzanova metoda, metoda bisekce nebo metoda dichotomie. Jedná se o nejjednodušší iterační metodu, která spočívá v utvoření posloupnosti intervalů Ik = [ak , bk ], Ik+1 ⊂ Ik , kde délka následujícího intervalu sk+1 je vždy poloviční oproti původnímu intervalu sk 1 1 sk+1 = sk = (bk − ak ), 2 2 tedy optimální bod x∗ ∈ Ik , a1 = a, b1 = b, k = 1 : N . Intuitivně je zřejmé, že lim Ik = x∗
k→∞
lim sk = 0.
k→∞
Z unimodálnosti funkce je zřejmé, že jestliže nalezneme hledané minimum x∗ na intervalu [a, b], pak nám tento bod rozdělí funkci f tak, že na intervalu [a, x∗ ] klesá a na intervalu [x∗ , b] roste. Z toho nám plyne, že ∀x ∈ [a, x∗ ] je f 0 (x) < 0 a ∀x ∈ [x∗ , b] je f 0 (x) > 0. Přirozeně platí, že f 0 (x∗ ) = 0. Na těchto vlastnostech postavíme metodu půlení intervalu.
V prvním kroku položíme 1 a1 = a, b1 = b, x1 = (a1 + b1 ), s1 = b1 − a1 = b − a, 2 a pro k = 1 : N opakujeme postup níže, dokud nedostaneme požadovanou přesnost 1. jestliže f 0 (xk ) < 0 ⇒ x∗ > xk , položíme ak+1 = xk , bk+1 = bk , 1 xk+1 = (ak+1 + bk+1 ) 2 12
2. jestliže f 0 (xk ) > 0 ⇒ x∗ < xk , položíme ak+1 = ak , bk+1 = xk , 1 xk+1 = (ak+1 + bk+1 ) 2 3. jestliže f 0 (xk ) = 0 ⇒ xk je hledané minimum. Lze odvodit, že nutný počet iterací pro dosažení požadované přesnosti ε je dán vztahem N ≥
log b−a 2ε . log 2
4.2. Metoda sečen Metoda sečen, někdy také nazývaná jako sekantová, vychází z Newtonovy metody. Počítání funkčních hodnot první a druhé derivace může být pracné, a tak aproximujeme f 00 (xk ) takto f 00 (xk ) ≈
f 0 (xk ) − f 0 (xk−1 ) . xk − xk−1
Následným dosazením do vzorce pro Newtonovu metodu dostaneme xk+1 = xk −
xk − xk−1 f 0 (xk ). 0 k ) − f (xk−1 )
f 0 (x
Takto jsme získali metodu sečen, která požaduje na vstup dvě počáteční iterace xk , xk−1 a funkční hodnoty první derivace v těchto bodech. Geometrický význam spočívá v nalezení iterace xk+1 jako průsečíku přímky [(xk , f (xk )), (xk−1 , f (xk−1 ))] s osou x.
4.3. Metoda regula falsi Metoda regula falsi, v literatuře někdy uváděná jako "method of false position", je dalším typem iterační metody, která vychází z metody Newtonovy. Newtonova metoda je založena na informacích o jediném bodu xk , funkční hodnotě a funkční hodnotě první a druhé derivace v tomto bodě. Užitím více bodů snížíme požadavky na každý z nich. Metoda regula falsi tedy spočívá v užití bodu xk a 13
v hodnotách f (xk ), f 0 (xk ), f 0 (xk−1 ). Poté můžeme vytvořit kvadratickou funkci h(x) h(x) = f (xk ) + f 0 (xk )(x − xk ) +
f 0 (xk−1 ) − f 0 (xk ) (x − xk )2 . , xk−1 − xk 2
která má odpovídající hodnoty stejné. Stejně jako u Newtonovy metody funkci zderivujeme, dosadíme do bodu xk+1 a položíme rovnu nule: h0 (x) = f 0 (xk ) +
h0 (xk+1 ) = f 0 (xk ) +
(f 0 (xk−1 ) − f 0 (xk ))(x − xk ) xk−1 − xk
(f 0 (xk−1 ) − f 0 (xk ))(xk+1 − xk ) =0 xk−1 − xk
A po úpravě dostáváme xk+1 = xk − f 0 (xk ).
xk−1 − xk 0 k−1 ) − f (xk )
f 0 (x
Metoda regula falsi také nezávisí na hodnotě f (xk ), proto naše aproximace povede buď přes f (xk ) nebo f (xk−1 ). Na tuto metodu se dá také nahlížet jako na aproximaci Newtonovy metody, kde druhou derivaci nahradíme rozdílem dvou prvních derivací ve dvou různých bodech. Další podobností s metodou tečen je, že se dá také psát ve tvaru f 0 (x) ≡ g(x) = 0, tedy také iterační předpis je možné přepsat na tvar xk+1 = xk − g(xk ).
14
xk−1 − xk . g(xk−1 ) − g(xk )
Obrázek 4: Princip metody regula falsi, první tři iterace
4.4. Metoda kubické interpolace Další metoda, která je při jednorozměrné minimalizaci užitečná, se nazývá metoda kubické interpolace. Je založena na polynomu třetího stupně neboli na kubickém polynomu: y = f (x) = a0 + a1 x + a2 x2 + a3 x3 . Z nákresu polynomu y je zřejmé, že polynom může mít minimum stejně tak jako maximum.
Obrázek 5: Kubická funkce, vlevo s koeficentem a3 > 0, vpravo s koeficientem a3 < 0 15
Nyní si ukážeme způsob, jakým nalezneme hledané minimum. Postupujeme jako při obvyklém hledání extrému, tedy položíme první derivaci funkce rovno nule: y 0 = f 0 (x) = a1 + 2a2 x + 3a3 x2 = 0 Vypočteme kořeny této kvadratické rovnice klasickým způsobem a dostáváme x=
−2a2 ±
p p 4a22 − 12a1 a3 −a2 ± a22 − 3a1 a3 = 6a3 3a3
Z matematické analýzy plyne, že v bodech lokálního minima je druhá derivace kladná, takže nám platí následující vztah pro druhou derivaci v hledaném bodě minima x∗ y 00 (x∗ ) = 2a2 + 6a3 x∗ > 0 neboli x∗ > −
a2 3a3
(1)
Za vztahu (1) a odvození je nyní zřejmé, který ze dvou kořenů je minimum. Neznámé koeficienty a0 , a1 , a2 , a3 se dají určit několika způsoby, jedním z nich je vycházet ze systému rovnic y1 = f (x1 ) = a0 + a1 x1 + a2 x21 + a3 x31 y2 = f (x2 ) = a0 + a1 x2 + a2 x22 + a3 x32 y3 = f (x3 ) = a0 + a1 x3 + a2 x23 + a3 x33 , a k tomuto systému přidat funkční hodnotu první derivace v jednom z bodů x1 , x2 , x3 , vybereme například bod x1 , tedy platí y 0 = f 0 (x) = a1 + 2a2 x + 3a3 x2 y10 = f 0 (x1 ) = a1 + 2a2 x1 + 3a3 x21 Řešením tohoto systému čtyř rovnic o čtyřech neznámých nalezneme hodnoty koeficientů a0 , a1 , a2 , a3 , přičemž koeficient a0 nás početně nezajímá. Hodnoty koeficientů jsou nicméně velmi rozsáhlé, proto je zapisujeme touto formou: 16
a1 = y10 − 2a2 x1 − 3a3 x21 a2 = α − γa3 α−β a3 = , γ−δ kde y2 − y1 + y10 (x1 − x2 ) (x1 − x2 )2 y3 − y1 + y10 (x1 − x3 ) β= (x1 − x3 )2 2 2x − x2 (x1 + x2 ) γ= 1 (x1 − x2 ) 2x21 − x3 (x1 + x3 ) δ= (x1 − x3 ) α=
Je zřejmé, že při užití metody kubické interpolace zahrnuje jedna iterace mnohem více počítání, než-li při užití interpolace kvadratické. I přesto však může být kubická interpolace účinnější než kvadratická, protože aproximace funkce f (x) polynomem třetího řádu je přesnější než aproximace polynomem řádu dva.
Obrázek 6: Příklad užití metody kubické interpolace, funkci f (x) = x5 − 2x proložíme kubickou křivkou y = 1.75x3 − 0.75x2 − 2x
17
5. Metody nevyužívající derivace 5.1. Metoda kvadratické interpolace Metoda kvadratické interpolace (někdy také známá jako Powellova metoda) je iterační metoda, která požaduje na vstup do algoritmu (a také do každé další iterace) tři body - (xi , yi ), yi = f (xi ), i = 1, 2, 3. Těmito body poté proložíme kvadratickou parabolu y = a0 + a1 x + a2 x 2
(2)
a najdeme její extrém, jehož souřadnice označíme (x4 , y4 ). Pro extrém kvadratické fce (2) platí dy a1 = a1 + 2a2 x => x4 = x = − . dx 2a2 Řešení systému rovnic y1 = a0 + a1 x1 + a2 x21 y2 = a0 + a1 x2 + a2 x22 y3 = a0 + a1 x3 + a2 x23 vede k hodnotám koeficientů a1 , a2 . Po odečtení první rovnice od druhé dostáváme výraz y2 − y1 = a1 (x2 − x1 ) + a2 (x22 − x21 ) => a1 =
y2 − y1 − a2 (x22 − x21 ) x2 − x 1
Po odečtení první rovnice od třetí dostáváme výraz y3 − y1 = a1 (x3 − x1 ) + a2 (x23 − x21 ) a po dosazení již odvozeného a1 dostaneme y3 − y1 =
y2 − y1 − a2 (x22 − x21 ) (x3 − x1 ) + a2 (x23 − x21 ) x2 − x1
(y3 − y1 )(x2 − x1 ) = (y2 − y1 − a2 (x22 − x21 ))(x3 − x1 ) + a2 (x23 − x21 )(x2 − x1 ) 18
(y3 − y1 )(x2 − x1 ) = (y2 − y1 )(x3 − x1 ) + a2 ((x23 − x21 )(x2 − x1 ) − (x22 − x21 )(x3 − x1 )) a2 =
(y3 − y1 )(x2 − x1 ) − (y2 − y1 )(x3 − x1 ) (x23 − x21 )(x2 − x1 ) − (x22 − x21 )(x3 − x1 )
Dosazením do rovnice x = x4 = −
a1 dostáváme výraz 2a2
−(y2 − y1 ) + a2 (x22 − x21 ) a1 x2 − x1 x=− = = 2a2 2a2
=
(x2 −x1 )(y3 −y1 )−(x3 −x1 )(y2 −y1 ) (x22 (x23 −x21 )(x2 −x1 )−(x22 −x21 )(x3 −x1 ) (x2 −x1 )(y3 −y1 )−(x3 −x1 )(y2 −y1 ) 2 (x 2 −x2 )(x −x )−(x2 −x2 )(x −x ) (x2 − x1 ) 2 1 3 1 3 1 2 1
−(y2 − y1 ) +
− x21 )
=
1 (x2 − x21 )(y3 − y1 ) − (x23 − x21 )(y2 − y1 ) = . 2 2 (x2 − x1 )(y3 − y1 ) − (x3 − x1 )(y2 − y1 ) Další postup při hledání extrému funkce f na intervalu [a, b]: 1. V bodech x1 = a, x2 = b, x3 =
a+b vypočítáme hodnoty funkce yi = 2
f (xi ), i = 1, 2, 3. 2. Využitím vztahu pro x4 vypočítáme další aproximaci polohy extrému. Dále vypočítáme y4 = f (x4 ). 3. Zjistíme, zda hodnoty funkce vyhovují požadavku na přesnost, tedy jestli platí vztah |x3 − x4 | ≤ ε, jestliže ano, přejdeme k následující bodu. Jestliže vztah není splněn, je třeba udělat tyto substituce x1 = x2 , y1 = y2 x2 = x3 , y2 = y3 x3 = x4 , y3 = y4 ,
a vrátit se k bodu 2. 19
4. Výpočet nyní můžeme ukončit, souřadnice extrému jsou (x4 , y4 ).
Obrázek 7: Princip metody kvadratické interpolace
Nevýhoda metody kvadratické interpolace je v tom, že někdy může selhat i při unimodální funkci, například tehdy, když není účelová funkce při hledání minima konvexní. Důvodem je to, že při prokládání bodů parabolou může vzniknout parabola konkávní a extrém, který nalezneme, nebude hledané minimum, ale naopak maximum funkce. Poté metoda kvadratické interpolace nepočítá správné hodnoty a je třeba užít jinou metodu. Příklad 1. Metodou kvadratické interpolace nalezněte minimum funkce f (x) = x5 − 2x na intervalu [−1, −0.5], kde přesnost ε = 0.01. Zaokrouhlujte na tři desetinná místa. V prvním kroku kvadratické interpolace volíme jako body x1 ,x2 krajní body a+b intervalu a jako bod x3 = . Poté vypočítáme hodnoty yi = f (xi ) pro i = 1 : 3 2 a hodnotu x4 . 20
x1 = −1, x2 = −0.5, x3 = −0.75, . x4 = −0.794, |x3 − x4 | = 0.044 > 0.01.
y1 y2 y3 y4
=1 = 0.969 = 1.263 . = 1.272
Přesnost nebyla splněna, je potřeba udělat dané substituce a opakovat postup. x1 = −0.5, x2 = −0.75, x3 = −0.794, . x4 = −0.804, |x3 − x4 | = 0.01 = 0.01.
y1 y2 y3 y4
= 0.969 = 1.263 = 1.272 . = 1.272
Nyní je splněna přesnost, ale je zřejmé, že výsledek není správný, protože funkční hodnota v bodě x4 je 1.272, což je více než funkční hodnota v krajním bodě intervalu x1 = −0.5, kde f (x1 ) = y1 = 0.969, tedy v bodě x4 určitě neleží minimum. Naopak se nám povedlo nalézt maximum dané funkce na tomto intervalu. Důvodem je to, že funkce není na našem intervalu konvexní, i přesto, že je unimodální. Z definice konvexity platí, že pro libovolnou dvojici bodů x1 , x2 ∈ I a libovolné λ ∈ [0, 1] platí vztah f (λx1 + (1 − λ)(x2 )) ≤ λf (x1 ) + (1 − λ)f (x2 ). Zvolíme-li tedy za x1 = −1, x2 = −0.5 a λ = 0.5, dostáváme 1 f (−0.75) ≤ (f (−1) + f (−0.5)) =⇒ 1.263 ≤ 0.985, 2 což ale určitě neplatí.
5.2. Rovnoměrná komparativní metoda Rovnoměrná komparativní metoda spočívá v rozdělení intervalu [a, b] pomocí bodů xk = a + k
b−a ,k = 1 : N N +1 21
a nalezení f (xm ) = min f (xk ). 1≤k≤N
Optimální bod x∗ se tedy nachází v intervalu [xm−1 , xm+1 ] a lze psát x∗ = x m ±
b−a N +1
U komparativních metod nazýváme určení hodnoty účelové funkce (ať už výpočtem či experimentálním měřením) experiment. Efektivnost metody je vyjádřena poměrem délky počátečního intervalu k délce posledního j-tého intervalu neurčitosti, tedy E=
b−a . lj
Pokud není předem dán počet experimentů N , lze z požadované přesnosti ε tato hodnota odvodit, neboť lj = xm+1 − xm−1 = 2.
⇒
b−a ≤ 2ε N +1
b−a ≤ε N +1
b − a ≤ ε(N + 1) b − a ≤ εN + ε b − a − ε ≤ εN ⇒N ≥
b−a −1 ε
22
Vztah pro efektivnost se dá také upravit následujícím způsobem E=
b−a b−a N +1 = 2(b−a) = . lj 2 N +1
Obrázek 8: Princip rovnoměrné komparativní metody
5.3. Metoda zlatého řezu Metoda zlatého řezu se řadí mezi postupné komparativní metody. Stejně jako u metody půlení intervalu vytvoříme posloupnost intervalů Ik = [ak , bk ], Ik+1 ⊂ Ik ⊂ . . . I2 ⊂ I1 = [a, b], kde délka intervalů s se neustále snižuje sk+1 < sk < · · · < s2 < s1 = b − a. Je tedy zřejmé, že opět lim Ik = x∗
k→∞
23
lim sk = 0.
k→∞
Metoda zlatého řezu je založena na rozdělení intervalu neurčitosti následujícím způsobem sk sk−1 = = r = konst. sk sk+1 a platí sk−1 = sk + sk+1 , k = 2 : N
Obrázek 9: Dělení intervalu metodou zlatého řezu
Dosazením do proměnné r dostaneme kvadratickou rovnici r2 − r − 1 = 0, která má dva reálné kořeny. Zajímá nás jen kladný kořen √ 1+ 5 . = 1, 618034, r= 2 který se nazývá zlatý řez. Platí pro něj následující vztahy 1 = r
√ √ 5−1 . 1 3− 5 . 1 1 = 0, 618034, 2 = = 0, 381966, 2 + = 1. 2 r 2 r r
Obrázek 10: k-tý krok dělení metodou zlatého řezu 24
V prvním kroce položíme a1 = a, b1 = b, x11 = a1 +
1 1 (b1 − a1 ), x21 = b1 − 2 (b1 − a1 ), s1 = b1 − a1 = b − a 2 r r
a pro k = 2 : N − 1 opakujeme tento postup 1. jestliže f (x1k−1 ) ≤ f (x2k−1 ), položíme ak = ak−1 , bk = x2k−1 , x1k = ak +
1 (bk − ak ), x2k = x1k−1 r2
2. jestliže f (x1k−1 ) > f (x2k−1 ), položíme ak = x1k−1 , bk = bk−1 , x1k = x2k−1 , x2k = bk −
1 (bk − ak ) r2
V posledním, N -tém kroce, užijeme tento postup 1. jestliže f (x1k−1 ) ≤ f (x2k−1 ), položíme aN = aN −1 , bN = x2N −1 , 1 xN = (aN + bN ) 2 2. jestliže f (x1k−1 ) > f (x2k−1 ), položíme aN = x1N −1 , bN = bN −1 , 1 xN = (aN + bN ). 2 Lze odvodit, že nutný počet iterací pro dosažení požadované přesnosti ε je dán log b−a 2ε vztahem N ≥ + 1. log r
5.4. Fibonacciho metoda Fibonacciho metoda, někdy též známá pod pojmem Kieferova, je modifikací metody zlatého řezu. Je také založena na posloupnosti intervalů Ik = [ak , bk ] a na práci s délkami těchto intervalů. Při snaze zkrátit délky intervalů neurčitosti využívá Fibonacciho číselné posloupnosti. Počítací technika je podobná jako v metodě zlatého řezu. Algoritmus postupu je součástí přílohy.
25
6. Příklady V této kapitole si pro ilustraci uvedeme jednoduché příklady na určité typy metod. Příklad 2. Newtonovou metodou, metodou půlení intervalu, metodou sečen a metodou kubické interpolace nalezněte minimum funkce f (x) = x3 − 5x + 3 na intervalu [1, 2], kde přesnost ε = 0.05. Zaokrouhlujte na čtyři desetinná místa. V metodách užívajících derivace je zapotřebí první a druhá derivace f (x), tedy f 0 (x) = 3x2 − 5, f 00 (x) = 6x.
Newtonova metoda Jako vhodné se jeví zvolit za počáteční iteraci některý bod intervalu, tedy například zvolíme x1 = 1. Dále spočítáme hodnoty f 0 (x1 ) = −2, f 00 (x1 ) = 6. Druhou iteraci spočítáme dosazením do vzorce Newtonovy metody a dostáváme x2 = x1 −
4 f 0 (x1 ) = . 00 f (x1 ) 3
Poté otestujeme přesnost |x2 − x1 | =
1 > 0.05. 3
Přesnost není splněna, je třeba pokračovat další iterací. 1 f 0 (x2 ) = , 3
f 00 (x2 ) = 8,
x3 = x2 −
f 0 (x2 ) 31 . = = 1.2917, 00 f (x2 ) 24
|x3 − x2 | = 0.0417 < 0.05.
Výpočet nyní můžeme ukončit, požadovaný výsledek je (x3 , f (x3 )) = (1.2917, −1.3033).
26
Metoda půlení intervalu Jako první je dobré spočítat si potřebný počet iterací N , který zaručuje dosažení přesnosti. 1 b−a log 2ε = 0.1 = 3.32 =⇒ N = 4. log2 log2
log N≥
V prvním kroce pokládáme a1 = a = 1, b1 = b = 2, x1 = 1.5. Spočítáme f 0 (x1 ) = 1.75 > 0 a algoritmem metody bisekce pokračujeme. a2 = a1 = 1, a3 = x2 = 1.25, a4 = a3 = 1.25,
b2 = x1 = 1.5, b3 = b2 = 1.5, b4 = x3 = 1.375,
x2 = 1.25, x3 = 1.375, x4 = 1.3125.
f 0 (x2 ) = −0.3125 < 0 . f 0 (x3 ) = 0.6719 > 0
Námi hledané minimum je tedy v bodě (x4 , f (x4 )) = (1.3125, −1.3015). Metoda sečen V metodě sečen je zapotřebí na vstup dvě iterace. Zvolíme proto krajní body našeho intervalu a dopočítáme hodnoty potřebné pro vzorec. x1 = 1, x2 = 2, f 0 (x1 ) = −2, f 0 (x2 ) = 7, x2 − x1 . .f 0 (x2 ) = 1.2222, 0 2 ) − f (x1 ) x3 − x2 . x4 = x3 − 0 .f 0 (x3 ) = 1.2759, f (x3 ) − f 0 (x2 ) x4 − x3 . x5 = x4 − 0 .f 0 (x4 ) = 1.2914, 0 f (x4 ) − f (x3 )
x3 = x2 −
f 0 (x
. |x3 − x2 | = 0.7778 > 0.05 . |x4 − x3 | = 0.0537 > 0.05 . |x4 − x3 | = 0.0155 < 0.05
V páté iteraci jsme nalezli výsledek odpovídající kritériím přesnosti, (x5 , f (x5 )) = (1.2914, −1.3033).
27
Metoda kubické interpolace Jelikož zadaná funkce je třetího stupně, není třeba počítat koeficienty v kubické interpolaci, stačí si pouze uvědomit, že a0 = 3, a1 = −5,a2 = 0,a3 = 1. Poté dosadíme do vzorce pro minimum
x=
a2 ±
p p (−3)(5) . a22 − 3a1 a3 =± = ±1.291. 3a3 3
Z teorie metody kubické interpolace ale víme, že musí platit následující vztah, ze kterého je jasné, který z kořenů je minimum. x∗ > −
a2 > 0 =⇒ x∗ = 1.291. 3a3
Příklad 3. Metodou kvadratické interpolace, rovnoměrnou komparativní metodou a metodou zlatého řezu nalezněte minimum funkce f (x) = 3x3 − 2x + 7 na intervalu [0, 1], kde přesnost ε = 0.1. Zaokrouhlujte na tři desetinná místa.
Metoda kvadratické interpolace V kvadratické interpolaci volíme jako body x1 ,x2 krajní body intervalu a jako a+b bod x3 = . Dále vypočítáme hodnoty yi = f (xi ) pro i = 1 : 3 a hodnotu x4 . 2 Pak otestujeme přesnost. x1 = 0, x2 = 1, x3 = 0.5, . x4 = 0.389, |x3 − x4 | = 0.111 > 0.1.
y1 y2 y3 y4
=7 =8 = 6.375 . = 6.399
Přesnost není splněna, je třeba udělat substituce a zopakovat daný postup.
28
x1 = 1, x2 = 0.5, x3 = 0.389, . x4 = 0.464, |x3 − x4 | = 0.075 < 0.1.
y1 y2 y3 y4
=8 = 6.375 = 6.399 . = 6.372
Splnili jsme požadavek přesnosti a tím nalezli minimum funkce, kterým je bod (x4 , f (x4 )) = (0.464, 6.372). Rovnoměrná komparativní metoda V této metodě je vhodné začít spočítáním počtu iterací N pro zaručení přesnosti. N≥
b−a − 1 = 9 =⇒ N = 9. ε
k(b − a) spočítáme devět iterací a jejich funkční hodnoty. N +1 Z funkčních hodnot pak vybereme tu nejmenší. Pomocí vzorce xk = a+
x1 x2 x3 x4 x5 x6 x7 x8 x9
= 0.1, = 0.2, = 0.3, = 0.4, = 0.5, = 0.6, = 0.7, = 0.8, = 0.9,
f (x1 ) = 6.803 f (x2 ) = 6.624 f (x3 ) = 6.481 f (x4 ) = 6.392 f (x5 ) = 6.375 f (x6 ) = 6.448 f (x7 ) = 6.629 f (x8 ) = 6.936 f (x9 ) = 7.387
Nejmenší funkční hodnota je f (x5 ) = 6.375 =⇒ minimum se nachází v bodě (x5 , f (x5 )) = (0.5, 6.375). Metoda zlatého řezu b−a . 2ε + 1 = 4.3445 =⇒ N = 5. V pěti krocích Opět začneme výpočtem N ≥ logr tedy aplikujeme metodu zlatého řezu podle popisu. log
29
1. a1 = a = 0, b1 = b = 1, x11 = 0.382, x21 = 0.618 f (x11 ) = 6.403 < f (x21 ) = 6.472 2. a2 = a1 = 0, b2 = x21 = 0.618, x12 = 0.236, x22 = 0.382 f (x12 ) = 6.567 > f (x21 ) = 6.403 3. a3 = x12 = 0.236, b3 = b2 = 0.618, x13 = 0.382, x23 = 0.472 f (x13 ) = 6.403 > f (x23 ) = 6.372 4. a4 = x13 = 0.382, b4 = b = 3 = 0.618, x14 = 0.472, x24 = 0.528 f (x14 ) = 6.372 < f (x24 ) = 6.386 5. a5 = a4 = 0.382, b5 = x24 = 0.528, x15 = x25 = 0.455 Minimum jsme našli v bodě (x5 , f (x5 )) = (0.455, 6.373). Následující příklady budou popsány velmi stručným, ilustračním způsobem. Příklad 4. Metodami využívajícími derivace nalezněte minimum funkce f (x) = 6x4 − 7sin(x) na intervalu [0, 1], kde přesnost ε = 0.01. Zaokrouhlujte na čtyři desetinná místa. f 0 (x) = 24x3 − 7cos(x), f 00 (x) = 72x2 + sin(x)
Newtonova metoda Za počáteční bod nelze zvolit nula, protože ve vzorci Newtonovy metody bychom dostali nulu ve jmenovateli. Jako počáteční bod tedy zvolíme x1 = 1. x2 x3 x4 x5
= 0.7404, = 0.6369, = 0.6167, = 0.6193,
|x2 − x1 | = 0.2596 > 0.01 |x3 − x2 | = 0.1035 > 0.01 |x4 − x3 | = 0.0172 > 0.01 |x5 − x4 | = 0.0004 < 0.01
Minimum nalezené Newtonovou metodou leží v bodě (0.6193, −3.1807). 30
Metoda půlení intervalu b−a 1 log 2ε = 0.02 = 5.64 =⇒ N = 6. log2 log2
log N≥
a1 a2 a3 a4 a5 a6
= a = 0, = x1 = 0.5, = a2 = 0.5, = a3 = 0.5, = x4 = 0.5625, = x5 = 0.5938,
b1 b2 b3 b4 b5 b6
= b = 1, = b1 = 1, = x2 = 0.75, = x3 = 0.625, = b4 = 0.625, = b5 = 0.625,
x1 x2 x3 x4 x5 x6
= 0.5 = 0.75, = 0.625, = 0.5625 = 0.5938, = 0.6094,
f 0 (x1 ) = −3.1431 < 0 f 0 (x2 ) = 5.0032 > 0 f 0 (x3 ) = 0.1826 > 0 f 0 (x4 ) = −1.645 < 0 f 0 (x5 ) = −0.7768 < 0
Námi hledané minimum je (0.6094, −3.1792).
Metoda sečen, metoda regula falsi Metoda regula falsi vychází v tomto příkladu při zaokrouhlení na čtyři desetinná místa stejně jako metoda sečen, proto jsou jejich hodnoty uvedeny společně. x1 = 0, x2 = 1, f 0 (x1 ) = −7, f 0 (x2 ) = 20.2179, x3 x4 x5 x6 x7 x8 x9
= 0.2572 = 0.4350 = 0.8260 = 0.5650 = 0.6049 = 0.6205 = 0.6193
|x3 − x2 | = 0.7428 > 0.01 |x4 − x3 | = 0.1958 > 0.01 |x5 − x4 | = 0.391 > 0.01 |x6 − x5 | = 0.261 > 0.01 |x7 − x6 | = 0.0399 > 0.01 |x8 − x7 | = 0.0156 > 0.01 |x9 − x8 | = 0.0012 < 0.01
31
Metoda kubické interpolace Za vstupní body zvolíme hodnoty x1 = 0,x2 = 0.5,x3 = 1. x1 = 0 x2 = 0.5 x3 = 1
y1 = 0 y2 = −2.981 y3 = 0.1097 y10 = −7 a3 = 10.0674 a2 = −2.9577 a1 = −7
α = 2.076 β = 7.1097 γ = 0.5 δ=1 xm1 = −0.3934 xm2 = 0.5892 −a2 > 0.0979 x> 3a3 x = xm2 = 0.5892 ym = max(y1 , y2 , y3 ) ⇒ m = 3 x0 = x = 0.5892 xm = x3 = x = 0.5892
x1 = 0 x2 = 0.5 x3 = 0.5892
y1 = 0 y2 = −2.981 y3 = −3.1668 y10 = −7 a3 = 7.6502 a2 = −1.7491 a1 = −7
α = 2.076 β = 2.7584 γ = 0.5 δ = 0.5892 xm1 = −0.4813 xm2 = 0.6337 −a2 > 0.0762 x> 3a3 x = xm2 = 0.6337 ym = max(y1 , y2 , y3 ) ⇒ m = 1 x0 = x = 0.6337 xm = x1 = x = 0.6337
32
x1 = 0.6337 x2 = 0.5 x3 = 0.5892
y1 = −3.1773 y2 = −2.981 y3 = −3.1668 y10 = 0.4666 a3 = 14.759 a2 = −11.6138 a1 = −2.5952
α = 14.4713 β = 15.7878 γ = 1.7674 δ = 1.8566 xm1 = −0.0947 xm2 = 0.6193 −a2 > 0.262 x> 3a3 x = xm2 = 0.6193
Aproximací kubickým polynomem jsme nalezli minimum (0.6193, −3.1807). Příklad 5. Metodami bez využití derivací nalezněte minimum funkce f (x) = sinh(x) − 4x2 na intervalu [−1, 0], kde přesnost ε = 0.1. Zaokrouhlujte na čtyři desetinná místa.
Metoda kvadratické interpolace x1 x2 x3 x4
= −1 =0 = −0.5 = −0.1217
y1 = 2.8248 y2 = 0 y3 = 0.4789 |x3 − x4 | = 0.3738 > 0.1
x1 x2 x3 x4
=0 = −0.5 = −0.1217 = −0.1271
y1 = 0 y2 = 0.4789 y3 = −0.0628 |x3 − x4 | = 0.0054 < 0.1
Minimum bylo nalezeno v bodě (−0.1271, −0.0628).
33
Rovnoměrná komparativní metoda N≥
x1 x2 x3 x4 x5 x6 x7 x8 x9
= −0.9, = −0.8, = −0.7, = −0.6, = −0.5, = −0.4, = −0.3, = −0.2, = −0.1,
b−a − 1 = 9 =⇒ N = 9. ε
f (x1 ) = 2.2135 f (x2 ) = 1.6719 f (x3 ) = 1.2014 f (x4 ) = 0.8033 f (x5 ) = 0.4789 f (x6 ) = 0.2292 f (x7 ) = 0.0555 f (x8 ) = −0.0413 f (x9 ) = −0.0602
Nejmenší funkční hodnota je f (x9 ), minimum se nachází v bodě (−0.1, −0.0602). Metoda zlatého řezu b−a . 2ε + 1 = 4.3445 =⇒ N = 5 logr
log N≥
1. a1 = −1, b1 = 0, x11 = −0.618, x21 = −0.382 f (x11 ) = 0.8696 > f (x21 ) = 0.1923 2. a2 = −0.618, b2 = 0, x12 = −0.382, x22 = −0.236 f (x12 ) = 0.1923 > f (x21 ) = −0.0154 3. a3 = −0.382, b3 = 0, x13 = −0.236, x23 = −0.1459 f (x13 ) = −0.0154 > f (x23 ) = −0.0613 4. a4 = −0.236, b4 = 0, x14 = −0.1459, x24 = −0.0902 f (x14 ) = −0.0613 < f (x24 ) = −0.0579 5. a5 = −0.236, b5 = −0.0902, x15 = x25 = −0.1631 34
Minimum jsme našli v bodě (x5 , f (x5 )) = (−0.1631, −0.0574). Fibonacciho metoda b−a √ . 2ε + log 5 − 1 = 4.016 =⇒ N = 5 logr logr
log N≥
F1 F2 F3 F4 F5
= 1.3416 =2 = 3.1305 =5 = 8.0498
1. a1 = −1, b1 = 0, x11 = −0.6111, x21 = −0.3889 f (x11 ) = 0.8439 > f (x21 ) = 0.2062 2. a2 = −0.6111, b2 = 0, x12 = −0.3889, x22 = −0.2444 f (x12 ) = 0.2082 > f (x21 ) = −0.0079 3. a3 = −0.3889, b3 = 0, x13 = −0.2444, x23 = −0.1667 f (x13 ) = −0.0079 > f (x23 ) = −0.0563 4. a4 = −0.2444, b4 = 0, x14 = −0.1667 = x24 f (x14 ) = −0.0563 > f (x14 + g) = −0.0566, kde g=0.001 5. a5 = −0.1667, b5 = 0, x15 = x25 = −0.0833 Minimum leží v bodě (−0.0833, −0.0556).
35
Literatura [1] Antoniou, A., Lu, W-S.: Practical Optimization: Alghoritms and Engineering applications. Springer, 2007. [2] Brunovská, A.: Malá optimalizácia: Metódy, programy, príklady. Alfa, Bratislava, 1990. [3] Luenberger, D.G., Ye, Y.: Linear and Nonlinear Programming. 3rd Edition, Springer, 2008. [4] Machalová, J.: Numerické metody (učební text), 2011. [5] Machalová, J., Netuka, H.: Numerické metody nepodmíněné optimalizace. PřF UP, 2013. [6] Míka, S.: Matematická optimalizace. Vydavatelství ZČU Plzeň, 1997. [7] Netuka, H.: Numerické metody optimalizace (učební text), 2012. [8] http://mi21.vsb.cz/sites/mi21.vsb.cz./files/unit/metody_optimalizace_obr.pdf [online 5. 1. 2013] (Dostál, Z., Beremlijski, P.: Metody optimalizace - interaktivní verze. VŠB - TUO, ZČU, 2012) [9] http://is.muni.cz/th/150978/prif_b/Novy3.pdf [online 5. 1. 2013] (Dokanová, L.: Konvexní funkce. PřF MU, 2007) [10] http://is.muni.cz/th/72828/prif_b/min.pdf [online 10. 12. 2012] (Urbánek, J.: Hledání minima funkce jedné proměnné. PřF MU, 2006) [11] http://cs.wikipedia.org/wiki/Metoda_tečen [online 10.2.2013] (Wikipedia.org: Metoda tečen) [12] http://books.fs.vsb.cz/StatickaOptimalizace [online 11.10.2012] (Vítečková, M., Jedlička, D.: Statická optimalizace systému. VŠB-TU Ostrava, 2003)
36
