ˇ MASARYKOVA UNIVERZITA V BRNE ˇ I´RODOVE ˇ DECKA ´ FAKULTA PR
ˇ ESˇENI´ NUMERICKE´ R ´ RNI´CH ROVNIC NELINEA ´ ZNI´ PRA ´ CE RIGORO
Mgr. Michal Sˇmerek Brno 2005
ii
Prohlášení: Prohlašuji, že předložená práce je mým původním autorským dílem, které jsem vypracoval samostatně. Veškerou literaturu a další zdroje, z nichž jsem při zpracování čerpal, v práci řádně cituji a jsou uvedeny v seznamu použité literatury.
V Brně dne .........................................
iii
ABSTRAKT Práce se zabývá studiem nelineárních rovnic, se zřetelnou specializací na polynomické rovnice. Autor předkládá řadu metod, jak určit kořeny polynomické rovnice. Podrobně studuje citlivost kořenů na malou změnu polynomu. Je všeobecně známo, jak získat kořeny polynomů, jak je získat s předepsanou přesností. Méně se ale klade důraz na přesné určení polynomu, resp. jeho koeficientů. Práce dokládá, že i malá nepřesnost, byť i jen v jednom koeficientu polynomu, může významným způsobem ovlivnit hodnoty kořenů. Text práce je rozdělen do tří částí. V první části jsou uvedeny základní metody řešení nelineárních rovnic, jako je metoda půlení, metoda regula falsi, metoda sečen, Newtonova metoda a rozšířená Newtonova metoda. U každé z nich je popsán algoritmus, jež určuje posloupnost aproximací, která konverguje ke hledanému řešení. Je definován řád konvergence metody, u jednotlivých metod je tento řád stanoven a většinou i dokázán. Srovnání jednotlivých metod z hlediska rychlosti konvergence je rovněž demonstrováno na testovacím příkladu. Druhá a třetí část práce se věnuje speciálnímu typu nelineárních rovnic, problematice polynomických rovnic. První z nich popisuje výpočet reálných kořenů polynomu libovolného stupně. Nejprve je vyložen postup, jak určit největší kořen daného polynomu, a následně je prezentován algoritmus pro získání všech ostatních kořenů polynomu. Je předložena dvojkroková metoda, která může nalezení kořenů polynomu značně urychlit ve srovnání s Newtonovou metodou. Práce obsahuje srovnání dvou metod pro určení všech kořenů polynomu, Maehlyho metody a metody dělení kořenovým činitelem. V posledním oddíle je pozornost zaměřena na citlivost kořenů polynomu. Je zde popsána závislost změny kořenů na malé změně některého z koeficientů polynomu. V souvislosti s tím je definován pojem špatně podmíněný polynom. Všechny části práce jsou proloženy názornými příklady. Vlastní výsledky autora jsou obsaženy ve druhé a třetí části práce. Hlubší studium této problematiky je nemyslitelné bez využití výpočetní techniky. Naopak, všechny v práci uváděné postupy jsou dobře algoritmizovatelné. Proto může studované téma dobře posloužit jako oblast matematiky vyučovaná s podporou počítače.
1
Obsah Předmluva a historický úvod 1 Metody řešení nelineárních rovnic 1.1 Metoda prosté iterace . . . . . . . 1.2 Metoda půlení . . . . . . . . . . . 1.3 Metoda regula falsi . . . . . . . . 1.4 Metoda sečen . . . . . . . . . . . 1.5 Metoda Newtonova . . . . . . . . 1.6 Rozšířená Newtonova metoda . . 1.7 Srovnání metod . . . . . . . . . .
2
. . . . . . .
5 6 9 12 18 20 23 26
2 Výpočet kořenů polynomu 2.1 Určení největšího kořene polynomu . . . . . . . . . . . . . . . . . 2.2 Určení ostatních kořenů polynomu . . . . . . . . . . . . . . . . . .
27 27 33
3 Citlivost kořenů polynomu 3.1 Špatně podmíněné polynomy . . . . . . . . . . . . . . . . . . . . . 3.2 Příklady polynomů druhého stupně . . . . . . . . . . . . . . . . . 3.3 Polynom vyššího stupně . . . . . . . . . . . . . . . . . . . . . . .
37 37 39 46
Závěr
50
Literatura
51
Přílohy Životopis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Publikační činnost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . CD-ROM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53 53 54 55
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
Předmluva a historický úvod
2
Předmluva a historický úvod Tématem práce je studium nelineárních rovnic, přičemž hlavní pozornost je věnována problematice řešení polynomických rovnic. Autor předkládá řadu metod, jak určit kořeny polynomické rovnice. Podrobně studuje citlivost kořenů na malou změnu polynomu. Je všeobecně známo, jak získat kořeny polynomů, jak je získat s předepsanou přesností. Méně se ale klade důraz na přesné určení polynomu, resp. jeho koeficientů. Práce dokládá, že i malá nepřesnost, byť i jen v jednom koeficientu polynomu, může významným způsobem ovlivnit hodnoty kořenů. Cílem práce bylo klasifikovat metody řešení nelineárních rovnic, popsat je a posoudit je z hlediska rychlosti konvergence. Těžištěm práce je ovšem problematika výpočtu kořenů polynomů a zkoumání citlivosti těchto kořenů na malé změny koeficientů polynomu. Text práce je rozdělen do tří částí. V první části jsou uvedeny základní metody řešení nelineárních rovnic, jako je metoda půlení, metoda regula falsi, metoda sečen, Newtonova metoda a rozšířená Newtonova metoda. U každé z nich je popsán algoritmus, jež určuje posloupnost aproximací, která konverguje ke hledanému řešení. Je definován řád konvergence metody, u jednotlivých metod je tento řád stanoven a většinou i dokázán. Srovnání jednotlivých metod z hlediska rychlosti konvergence je rovněž demonstrováno na testovacím příkladu. Druhá a třetí část práce se věnuje speciálnímu typu nelineárních rovnic, problematice polynomických rovnic. První z nich popisuje výpočet reálných kořenů polynomu libovolného stupně. Nejprve je vyložen postup, jak určit největší kořen daného polynomu, a následně je prezentován algoritmus pro získání všech ostatních kořenů polynomu. Je předložena dvojkroková metoda, která může nalezení kořenů polynomu značně urychlit ve srovnání s Newtonovou metodou. Práce obsahuje srovnání dvou metod pro určení všech kořenů polynomu, Maehlyho metody a metody dělení kořenovým činitelem. V poslední kapitole je pozornost zaměřena na citlivost kořenů polynomu. Je zde popsána závislost změny kořenů na malé změně některého z koeficientů polynomu. V souvislosti s tím je definován pojem špatně podmíněný polynom. V tomto ohledu je zvlášť diskutován případ jednoduchého kořene a násobného kořene. Ukázalo se totiž, že při velmi malé změně koeficientu polynomu probíhá změna jednotkového kořene kvalitativně odlišně než změna násobného kořene. Všechny části práce jsou proloženy názornými příklady. Práce obsahuje 26 obrázků a 21 tabulek. Byly zpracovány systémem MATLAB. Vlastní výsledky autora jsou obsaženy ve druhé a třetí části práce. Hlubší studium této problematiky je nemyslitelné bez využití výpočetní techniky. Naopak, všechny v práci uváděné postupy jsou dobře algoritmizovatelné. Dokladem toho je množství souborů na nosiči CD-ROM, jež je nedílnou součástí předkládané práce, viz příloha. Na CD-ROMu lze, kromě elektronické verze samotné práce, nalézt i programy, jež obsahují většinu v práci použitých algoritmů (m-soubory). Dále CD-ROM obsahuje také grafické výstupy (fig-soubory). Studované téma proto může dobře posloužit jako oblast matematiky vyučovaná
Předmluva a historický úvod
3
s podporou počítače. Zajímavý je i pohled na historický vývoj tohoto oboru (podrobněji viz článek [9]). Problém řešení polynomické rovnice p(x) = 0,
(0.1)
p(x) = an xn + an−1 xn−1 + · · · + a1 x + a0 ,
(0.2)
kde je polynom s reálnými koeficienty, byl znám již Sumerům ve 3. tis. před n. l. a hluboce ovlivnil rozvoj matematiky během mnoha staletí. Studium tohoto problému mělo velký význam pro rozvoj abstraktního myšlení a matematického zápisu. Motivovalo vznik a rozvoj mnoha odvětví matematiky, jako je teorie komplexních čísel, algebra, numerická matematika. Hlavním krokem v historii studia polynomických rovnic bylo uvážení problému v obecném abstraktním tvaru (0.2). Tento krok trval tisíciletí a vedl k zavedení moderního matematického formalismu. Mezitím, počínaje Sumerským a Babylonským obdobím, se studium zaměřilo na rovnice nižších řádů se speciálními koeficienty. Řešení speciálních kvadratických rovnic Babyloňany (kolem 2000 př. n. l.) a Egypťany, nalezené ve Rhindově a Ahmesově papyru (z 2. tis. př. n. l.), odpovídá užití našeho školního vzorce p −a1 ± a21 − 4a2 a0 . (0.3) x1,2 = 2a2 Pro plné pochopení tohoto vzorce je třeba znát teorii záporných, iracionálních a komplexních čísel. Důležitý úspěch v této oblasti byl formální důsledný důkaz Pythagorejců (kolem 500 př. n. l. ve starověkém Řecku), že rovnice x2 = 2 nemá racionální řešení, že řešení musí užít odmocnin a ne jen aritmetické operace. Snahy nalézt formuli, která by, podobně jako (0.3), zahrnovala jen aritmetické operace a odmocniny, uspěly v 16. století pro polynomy stupňů 3 a 4 (Scipione del Ferro, Nicolo Tartaglia, Ludovico Ferrari, Geronimo Cardano), ale hluboký vliv na matematiky mělo především selhání všech pokusů nalézt takovou formuli pro jakýkoli polynom stupně většího než 4. Přesněji řečeno takové pokusy vedly k větě, získané Ruffinim v roce 1813 a Abelem 1827, o neexistenci takové formule pro třídu polynomů stupně n pro žádné n > 4 a ke Galoisově fundamentální teorii z r. 1832. Tato teorie obsahuje důkaz neexistence řešení rovnice (0.2) ve tvaru formule s odmocninami již pro jednoduché speciální polynomické rovnice s přirozenými koeficienty, takové jako x5 − 4x − 2 = 0. Navzdory absenci řešící formule obecně platí základní věta algebry: Rovnice (0.1) má vždy (komplexní) řešení pro každý polynom p(x) každého kladného stupně n. Různé verze této věty podali Roth (1608), Girard (1629) a Descartes (1637), dále potom Rahn (1659), Newton (1685) a Maclaurin. Dokázat se povedla až v 19. století. Několik neúplných či chybných důkazů bylo již předtím
Předmluva a historický úvod
4
podáno D’Alembertem, Eulerem, Lagrangem a Gaussem (v jeho doktorské disertační práci v r. 1799). Chybu v Gaussově důkazu odstranil Ostrowski v r. 1920. Je snadné rozšířit základní větu algebry prokázáním existence faktorizace n Q p(x) = an (x − ξi ) pro každý polynom p(x), kde an 6= 0, takový, že ξ1 , ξ2 , . . . , ξn i=1
je n kořenů polynomu p(x) (ne nutně všechny rozdílné) a jsou jediným řešením rovnice (0.1). Po zjištění, že neexistuje přesná řešící formule, úsilí matematiků přešlo do navrhování iteračních algoritmů pro aproximaci řešení. V současné době seznam iteračních algoritmů navržených pro aproximaci řešení ξ1 , ξ2 , . . . , ξn rovnice (0.1) obsahuje stovky položek zahrnující období 4 tisíc let. Ve skutečnosti není důležité, zda počítač řešení získá pomocí formule nebo ne. Řešení je totiž zpravidla iracionální a tedy stejně nemůže být počítáno úplně přesně. S rozvojem vědy a techniky, se vznikem mnoha nových aplikací, vyvstává potřeba řešit polynomické rovnice řádu několika stovek či tisíců. Navíc se často požaduje vysoká přesnost stovek koeficientů i vlastního řešení. V těchto případech řešení rovnice (0.1) je problematické, a to motivuje další výzkum nových efektivních algoritmů pro řešení rovnice (0.1). Úsilí a zájem vědců v tomto směru je takový, že nové algoritmy pro řešení (0.1) se stále objevují každý rok.
Metody řešení nelineárních rovnic
1
5
Metody řešení nelineárních rovnic
Hledejme reálné kořeny rovnice f (x) = 0,
(1.1)
kde x je reálná proměnná a f (x) je v nějakém smyslu „rozumnáÿ funkce. Především budeme požadovat, aby funkce f (x) byla spojitá na nějakém intervalu. Předpokládejme, že kořen ξ rovnice (1.1) nelze nalézt v algebraickém tvaru. Musíme tedy hledat metody, které vedou k přibližnému řešení. Při řešení rovnice (1.1) můžeme užít následující věty. Věta 1.1 Nechť f je funkce spojitá na intervalu ha, bi a nechť v koncových bodech tohoto intervalu nabývá hodnot s opačnými znaménky, tj. f (a)f (b) < 0. Pak v intervalu (a, b) leží aspoň jeden kořen rovnice (1.1). Jestliže navíc existuje derivace f 0 funkce f a má konstantní znaménko v intervalu (a, b), pak funkce f (x) má v intervalu (a, b) právě jeden kořen ξ. Důkaz: Funkce f je spojitá na intervalu ha, bi. Nechť bez újmy na obecnosti je např. f (a) < 0 < f (b). Potom platí ∀h ∈ hf (a), f (b)i ∃x ∈ ha, bi : f (x) = h. Speciálně 0 ∈ (f (a), f (b)) ⇒ ∃ξ ∈ (a, b) : f (ξ) = 0. Tím je dokázána první část věty. Dále předpokládáme, že je buď f 0 (x) > 0 ∀x ∈ (a, b), anebo f 0 (x) < 0 ∀x ∈ (a, b). Nechť ξ1 , ξ2 ∈ (a, b) jsou kořeny funkce f , tzn. f (ξ1 ) = 0, f (ξ2 ) = 0. Pak je buď ξ1 = ξ2 , anebo ξ1 6= ξ2 . Pokud by bylo ξ1 6= ξ2 , pak z věty o střední hodnotě platí, že ∃x ∈ (ξ1 , ξ2 ) : f 0 (x) = 0, což je spor s předpokladem. Celkově, je-li f spojitá na ha, bi, f 0 (x) > 0 (resp. f 0 (x) < 0) ∀x ∈ (a, b), a platí-li f (a)f (b) < 0, pak v intervalu ha, bi existuje právě jeden kořen. 2 Všechny metody uvedené v této kapitole budeme aplikovat na příkladu nalezení nejmenšího kladného kořene funkce f (x) = cos x +
x − 1. 2
Graf této nelineární funkce je znázorněn na obrázku 1.
(1.2)
Metody řešení nelineárních rovnic
6
Obr. 1: Graf funkce f (x) = cos(x) + x/2 − 1
1.1
Metoda prosté iterace
Nechť f (x) je spojitá reálná funkce, která má tolik derivací, kolik jich bude v dalším třeba. Problém nalezení kořene rovnice (1.1) můžeme převést na ekvivalentní problém nalezení pevného bodu funkce F , tj. řešení rovnice x = F (x).
(1.3)
Pevný bod budeme označovat ξ. Věta 1.2 Nechť F je taková spojitá funkce, že ∀x ∈ ha, bi : F (x) ∈ ha, bi. Pak funkce F má v intervalu ha, bi pevný bod. Jestliže funkce F navíc splňuje Lipschitzovu podmínku, tj. |F (x) − F (y)| ≤ q |x − y| , ∀x, y ∈ ha, bi, 0 ≤ q < 1, pak F má v intervalu ha, bi jediný pevný bod. Důkaz: Jestliže F (a) = a, je a pevným bodem funkce F . Jestliže F (b) = b, je b pevným bodem funkce F . Předpokládejme dále, že F (a) > a, F (b) < b. Uvažujme funkci g, g(x) = F (x) − x. Zřejmě g(x) je spojitá funkce na intervalu ha, bi a platí g(a) = F (a) − a > 0, g(b) = F (b) − b < 0.
Metody řešení nelineárních rovnic
7
Z věty 1.1 plyne, že existuje bod ξ ∈ (a, b) : g(ξ) = 0, tj. F (ξ)−ξ = 0 ⇒ F (ξ) = ξ, a tedy ξ je pevným bodem funkce F (x). Nechť F (x) splňuje Lipschitzovu podmínku s konstantou q, 0 ≤ q < 1. Předpokládejme, že existují dva různé pevné body ξ, η ∈ ha, bi. Platí |ξ − η| = |F (ξ) − F (η)| ≤ q|ξ − η| < |ξ − η|, 2
což je spor a tedy ξ = η.
Důsledek 1.3 Nechť F je taková spojitá funkce, že ∀x ∈ ha, bi : F (x) ∈ ha, bi. Nechť v každém bodě x ∈ ha, bi existuje derivace F 0 (x) splňující podmínku |F 0 (x)| ≤ q < 1, pak F má v intervalu ha, bi jediný pevný bod. Důkaz: Plyne z věty o střední hodnotě: |F (x) − F (y)| = |F 0 (α)| |x − y| ≤ q |x − y| , α ∈ ha, bi. 2 Nechť jsou splněny předpoklady věty 1.2. Zvolme libovolnou počáteční aproximaci x1 ∈ ha, bi. Generujme posloupnost {xi }∞ i=1 následovně: xi+1 = F (xi ),
i = 1, 2, 3, . . .
(1.4)
Funkce F se nazývá iterační funkce a metoda (1.4) se nazývá iterační metoda nebo také metoda prosté iterace. Uvedená iterační metoda (1.4) je metodou jednobodovou, neboť výpočet aproximace xi+1 závisí pouze na předchozí aproximaci xi . Obecnější jsou potom n-bodové iterační metody, které mají tvar xi+1 = F (xi , xi−1 , . . . , xi−n+1 ),
i = n, n + 1, n + 2, . . . .
(1.5)
Aproximace xi+1 , daná vztahem (1.5), vždy závisí na n bezprostředně předchozích aproximacích. V tom případě hovoříme o stacionární iterační metodě. Ještě obecněji lze uvažovat nestacionární n-bodovou iterační metodu, xi = F (xk1 , xk2 , . . . , xkn ),
i = n + 1, n + 2, . . . ,
i > k1 > k2 > · · · > kn . (1.6)
Zde aproximace xi závisí vždy na n předchozích aproximacích, avšak ne nutně bezprostředně předcházejících aproximaci xi . Jak uvidíme dále, metoda půlení (viz odstavec 1.2) či metoda regula falsi (viz odstavec 1.3) jsou příklady dvoubodových nestacionárních metod; metoda sečen (viz odstavec 1.4) je dvoubodovou stacionární metodou; metoda Newtonova (viz odstavec 1.5) a rozšířená Newtonova metoda (viz odstavec 1.6) jsou zástupci jednobodových (stacionárních) metod. Otázkou nyní je, za jakých předpokladů a jak rychle bude posloupnost (1.4), resp. (1.5), resp. (1.6) konvergovat k pevnému bodu ξ. Právě z hlediska rychlosti konvergence má význam zavedení pojmu řád iterační metody.
Metody řešení nelineárních rovnic
8
Definice 1.4 Předpokládejme, že iterační metoda (1.4), resp. (1.5), resp. (1.6) je konvergentní, tj. lim xi = ξ. i→∞ Existuje-li reálné číslo p ≥ 1 takové, že platí |xi+1 − ξ| p = C 6≡ 0, i→∞ |xi − ξ| lim
říkáme, že iterační metoda je řádu p. Konstanta C se nazývá asymptotickou konstantou chyby a závisí na iterační funkci F . Požadavek C 6≡ 0 značí, že C 6= 0 pro obecnou funkci F . Tento požadavek zaručuje jednoznačnost čísla p. Je-li p = 1, musí být C ≤ 1, aby metoda konvergovala. V tom případě říkáme, že metoda konverguje lineárně. Jestliže pro nějakou funkci F je C = 0, pak daná iterační metoda konverguje rychleji než obvykle. Věta 1.5 Nechť funkce F má v okolí bodu ξ derivace až do řádu p ≥ 1 včetně. Iterační metoda (1.4) je řádu p, právě tehdy, když ξ = F (ξ),
F (j) (ξ) = 0, j = 1, 2, . . . , p − 1 ,
F (p) (ξ) 6= 0.
Důkaz: Taylorův rozvoj funkce F (x) v okolí bodu ξ je F (xi ) = ξ + (xi − ξ)F 0 (ξ) + . . . + = ξ+
(xi − ξ)(p−1) (p−1) (xi − ξ)p (p) F (ξ) + F (α) = (p − 1)! p!
(xi − ξ)p (p) F (α), p!
(1.7)
kde bod α leží v intervalu určeném body xi a ξ. Protože xi+1 = F (xi ), dostáváme ze vztahu (1.7) (xi − ξ)p (p) xi+1 − ξ = F (α), (1.8) p! a odtud
(p) F (ξ) |xi+1 − ξ| lim = 6≡ 0 p i→∞ |xi − ξ| p!
Metoda je tedy řádu p, p ∈ N. Naopak, nechť ∃j, 1 ≤ j < p : F (j) (ξ) 6= 0. Pak z (1.7) plyne, že metoda nemůže být řádu p. Stejně tak, jestliže platí F (p) (ξ) = 0, pak z (1.8) plyne, že metoda není řádu p. 2 Věta 1.6 Nechť F je taková spojitá funkce, že ∀x ∈ ha, bi : F (x) ∈ ha, bi. Nechť funkce F splňuje Lipschitzovu podmínku, tj. |F (x) − F (y)| ≤ q |x − y| , ∀x, y ∈ ha, bi, 0 ≤ q < 1. Pak pro libovolnou počáteční aproximaci x1 ∈ ha, bi je posloupnost {xi }∞ i=1 , xi = F (xi−1 ), konvergentní a platí lim xi = ξ, kde ξ je pevný bod funkce F . i→∞
Metody řešení nelineárních rovnic
9
Důkaz: Platí F (x) ∈ ha, bi ∀x ∈ ha, bi, a tedy pro x1 ∈ ha, bi je posloupnost {xi }∞ i=1 , xi = F (xi−1 ) definovaná a xi ∈ ha, bi pro i = 1, 2, . . . . Dále platí |xi − ξ| = |F (xi−1 ) − F (ξ)| ≤ q |xi−1 − ξ| . Odtud indukcí dostaneme |xi − ξ| ≤ q i−1 |x1 − ξ| . Protože 0 ≤ q < 1, platí lim |xi − ξ| = 0,
i→∞
a tedy posloupnost {xi }∞ i=1 konverguje k pevnému bodu ξ.
2
V následujících odstavcích se s některými iteračními metodami seznámíme blíže.
1.2
Metoda půlení
Metoda půlení neboli metoda bisekce je jednoduchou numerickou metodou pro výpočet kořene rovnice (1.1). Je založena na tvrzení věty 1.1. Metodu nyní popíšeme. Nechť jsou splněny předpoklady věty 1.1, tj. f je spojitá funkce na intervalu ha, bi. Předpokládejme dále, že na intervalu ha, bi má funkce f právě jeden kořen. Položíme a1 = a, b1 = b, x1 = 21 (a1 + b1 ). Nastane jedna ze tří možností: 1. Je-li f (a1 )f (x1 ) < 0, leží kořen ξ v intervalu ha1 , x1 i, položíme a2 = a1 , b2 = x1 . 2. Je-li f (x1 )f (b1 ) < 0, leží kořen ξ v intervalu hx1 , b1 i, položíme a2 = x1 , b2 = b1 . 3. Je-li f (x1 ) = 0, je x1 = ξ a kořen je nalezen. Položíme a2 = a3 = · · · = ξ a b2 = b3 = · · · = ξ. V prvních dvou případech jsou splněny předpoklady věty 1.1, opakujeme postup pro interval ha2 , b2 i. V každém případě získáme posloupnost intervalů ha1 , b1 i ⊃ ha2 , b2 i ⊇ . . . hai , bi i ⊇ . . . , přičemž f (ai )f (bi ) ≤ 0, i = 1, 2, . . . . Platí: a1 ≤ a2 ≤ · · · ≤ ξ, b1 ≥ b2 ≥ · · · ≥ ξ. Délky intervalů jsou dány vztahem b i − ai =
b−a , i = 1, 2, . . . , 2i−1
(1.9)
Metody řešení nelineárních rovnic
10
případně bi − ai = 0, i = k, k + 1, . . . ,
(1.10)
v případě, že jsme v k-tém kroku našli kořen ξ, tj. nastala třetí možnost. Protože ∀i ∈ N : a ≤ ai ≤ ξ ≤ bi ≤ b a lim (bi − ai ) = 0, platí i→∞
lim ai = lim bi = ξ.
i→∞
i→∞
Ukážeme, že ξ je kořenem rovnice f (x) = 0. Funkce f je spojitá a platí f (ai )f (bi ) ≤ 0, i = 1, 2, . . . . Dále lim f (ai )f (bi ) = lim f (ai ) lim f (bi ) = f 2 (ξ) ≤ 0.
i→∞
i→∞
i→∞
Odtud plyne, že f (ξ) = 0. Věta 1.7 Metoda půlení generuje posloupnost {xi }∞ i=1 , xi = guje ke kořenu ξ, přičemž chyba i-té aproximace je |xi − ξ| ≤
b−a , i = 1, 2, . . . . 2i
ai +bi , 2
která konver-
(1.11)
Důkaz: Platnost věty plyne z předchozích úvah a z faktu, že kořen ξ leží v intervalu hai , bi i, jehož délka je dána vztahem (1.9), příp. (1.10), a ze skutečnosti, že aproximace xi půlí tento interval. 2 Poznámka 1.8 Metodu půlení lze ekvivalentním způsobem zapsat také pomocí vztahu xi−1 + xs xi = , i = 1, 2, . . . , (1.12) 2 kde označíme x−1 = a, x0 = b a s = s(i), −1 ≤ s ≤ i − 2 je největší index takový, že f (xi−1 )f (xs ) < 0. Poznámka 1.9 Metoda půlení má řád konvergence roven jedné. Pokud vyloučíme singulární možnost, že v k-tém kroku najdeme kořen ξ, můžeme odhadnout , viz vztah (1.11), tj. (shora) chybu aproximace xi výrazem b−a 2i |xi − ξ| ≈ Potom
b−a , i = 1, 2, . . . . 2i
|xi+1 − ξ| 1 ≈ , |xi − ξ| 2
(1.13)
což ukazuje na to, že řád metody půlení je roven jedné. Příklad 1.1 Metodou půlení určete kořen funkce f (x) = cos(x) + x/2 − 1 ležící v intervalu h π4 , π2 i.
Metody řešení nelineárních rovnic
11
Snadno lze určit hodnoty funkce f v krajních bodech daného intervalu, tj. . . f (π/4) = 0,100; f (π/2) = −0,215. Rovněž lze snadno vyšetřit, že derivace f 0 (x) = − sin x + 1/2 nabývá v intervalu h π4 , π2 i pouze záporných hodnot. Z věty 1.1 tedy plyne, že v daném intervalu skutečně existuje právě jeden kořen funkce f . Vlastní výpočet metodou půlení je proveden v tabulce 1 postupem popsaným výše. Graficky je tento proces znázorněn na obrázku 2.
Obr. 2: Metoda půlení, f (x) = cos(x) + x/2 − 1, a = π/4, b = π/2
i
ai
bi
1
1 π 4 1 π 4 5 π 16 11 π 32 11 π 32 45 π 128 45 π 128 45 π 128
1 π 2 3 π 8 3 π 8 3 π 8 23 π 64 23 π 64 91 π 256 181 π 512
2 3 4 5 6 7 8
xi 3 π 8 5 π 16 11 π 32 23 π 64 45 π 128 91 π 256 181 π 512 361 π 1024
= 1,178 097 = 0,981 748 = 1,079 922 = 1,129 010 = 1,104 466 = 1,116 738 = 1,110 602 = 1,107 534
b i − ai 1 π 4 1 π 8 1 π 16 1 π 32 1 π 64 1 π 128 1 π 256 1 π 512
Tab. 1: Metoda půlení, f (x) = cos(x) + x/2 − 1, a = π/4, b = π/2
Metody řešení nelineárních rovnic
12
. 361 Z tabulky plyne, že x8 = 1024 π = 1,107 534. . 1 π = 0,003 068. Hodnota funkce v tomto bodě Chyba aproximace je |x8 − ξ| ≤ 1024 . je f (x8 ) = 0,000 636.
1.3
Metoda regula falsi
Předpokládejme i zde, že f je spojitá funkce na intervalu ha, bi, pro níž platí f (a)f (b) < 0. Metoda půlení jako další aproximaci užívá středu intervalu ha, bi, bez ohledu na to, jakých hodnot nabývá funkce f v krajních bodech tohoto intervalu. Jiná situace je u metody regula falsi - užívá nulového bodu přímky vedené body [a, f (a)], [b, f (b)], viz obrázek 3.
Obr. 3: Metoda regula falsi
Z obrázku je zřejmé (na základě podobnosti trojúhelníků), že platí b−x b−a = , f (b) f (b) − f (a)
(1.14)
kde bod x označuje nulový bod dané přímky. Ze vztahu (1.14) snadno vyjádříme hodnotu x b−a x=b− f (b). (1.15) f (b) − f (a) Položíme-li a = x1 , b = x2 , můžeme pomocí vztahu (1.15) vyjádřit následující
Metody řešení nelineárních rovnic
13
iteraci
x2 − x1 f (x2 ). (1.16) f (x2 ) − f (x1 ) Pro další výpočet je potřeba vzít v úvahu bod x3 a ten z bodů xi (i = 1, 2), který splňuje podmínku f (x3 )f (xi ) < 0. Obecně můžeme iteraci xi+1 počítat podle vzorce xi − xs xi+1 = xi − f (xi ), i = 2, 3, . . . , (1.17) f (xi ) − f (xs ) x3 = x2 −
kde s = s(i) < i je největší index, pro něž platí f (xi )f (xs ) < 0. Věta 1.10 Nechť f je funkce spojitá na intervalu ha, bi a nechť v koncových bodech tohoto intervalu nabývá hodnot s opačnými znaménky, tj. f (a)f (b) < 0. Nechť v intervalu (a, b) leží právě jeden kořen ξ. Nechť funkce f má na intervalu ha, bi spojité derivace f 0 (x) a f 00 (x). Dále nechť derivace f 0 (x), ani f 00 (x) nemění znaménko v intervalu ha, bi. Pak metoda regula falsi generuje posloupnost {xi }∞ i=1 , která konverguje ke kořenu ξ. Řád konvergence metody regula falsi je roven jedné. Důkaz: Předpokládejme např.: 1) f 0 (x) < 0 ∀x ∈ ha, bi, 2) f 00 (x) > 0 ∀x ∈ ha, bi, viz situace znázorněná na obrázku 4.
Obr. 4: Metoda regula falsi
Metody řešení nelineárních rovnic
14
Předpoklad 1) spolu s podmínkou f (a)f (b) < 0 implikuje f (a) > 0, f (b) < 0. Zvolme počáteční aproximace x1 = a, x2 = b. Bod x3 je definován vztahem (1.16). Platí pro něj x2 > x3 > ξ. Přitom, první nerovnost plyne přímo ze vztahu (1.16) uvědomíme-li si, že x2 − x1 > 0, f (x2 ) − f (x1 ) < 0, f (x2 ) < 0; druhá nerovnost plyne z toho, že f je funkce konvexní a klesající, a tedy průsečík sečny s osou x leží vpravo od kořene ξ. Dále, f (x3 ) < 0, takže indukcí plyne obecný vztah xi+1 = xi −
xi − x1 f (xi ), f (xi )f (x1 ) < 0, i = 2, 3, . . . . f (xi ) − f (x1 )
(1.18)
To znamená, že všechny sečny prochází bodem [x1 , f (x1 )]. Posloupnost {xi }∞ i=2 zřejmě monotónně konverguje k bodu ξ: x2 > x3 > · · · > xi > xi+1 > · · · > ξ, lim xi = ξ. i→∞
Pro tento případ ukážeme, že řád metody je roven jedné. Odečteme-li od obou stran rovnice (1.18) hodnotu ξ, dostáváme xi+1 − ξ = xi − ξ −
xi − x1 f (xi ), i = 2, 3, . . . . f (xi ) − f (x1 )
Odtud po úpravě xi+1 − ξ xi − x1 f (xi ) 1 = 1− = 1− k(xi , ξ), i = 2, 3, . . . , (1.19) xi − ξ f (xi ) − f (x1 ) xi − ξ k(xi , x1 ) kde k(xi , x1 ) je směrnice sečny si1 procházející body [xi , f (xi )], [x1 , f (x1 )]; k(xi , ξ) je směrnice sečny siξ procházející body [xi , f (xi )], [ξ, 0]. Jelikož posloupnost {xi }∞ i=2 konverguje k bodu ξ, tj. lim xi = ξ, potom, v limitním i→∞
přechodě i → ∞, přejde sečna si1 v sečnu s procházející body [ξ, 0] a [x1 , f (x1 )], sečna siξ přejde v tečnu t funkce f v bodě ξ, viz obrázek 4. Můžeme tedy vztah (1.19) přepsat takto: xi+1 − ξ 1 =1− f 0 (ξ) = c, i→∞ xi − ξ k(ξ, x1 ) lim
(1.20)
kde c je konstanta, c ∈ (0, 1), viz poznámka 1.11. Podobně bychom poslední vztah odvodili i pro ostatní případy. Z uvedeného plyne, že metoda regula falsi má za daných předpokladů řád konvergence roven jedné. 2 Poznámka 1.11 Ukážeme, že c ve vztahu (1.20) je konstantou z intervalu (0, 1). Směrnice k(ξ, x1 ) sečny s a směrnice f 0 (ξ) tečny t mají stejná znaménka, neboť platí předpoklad 1), tj. f 0 (x) < 0 ∀x ∈ ha, bi, a body x1 i ξ leží v intervalu ha, bi. Z věty o střední hodnotě plyne, že ∃η ∈ (x1 , ξ) : f 0 (η) = k(ξ, x1 ).
Metody řešení nelineárních rovnic
15
Jelikož a < η < ξ < b, a z platnosti předpokladů 1) a 2) plyne |k(ξ, x1 )| = |f 0 (η)| > |f 0 (ξ)|. (Tzn., že sečna s je vždy strmější než tečna t.) 1 Odtud už je zřejmé, že c = 1 − f 0 (ξ) = konst. ∈ (0, 1). k(ξ, x1 ) Na tomto místě ještě poznamenejme, že konvergence je tím rychlejší, čím je konstanta c blíže číslu nula. Pokud je číslo c blízké číslu jedna, je naopak konvergence velice pomalá. Tuto situaci demonstruje příklad 1.3. Poznámka 1.12 Předpoklad o konstantním znaménku funkcí f 0 (x) a f 00 (x) na intervalu ha, bi je ve větě nadbytečný. Důkaz příslušné obecnější (bez tohoto předpokladu) věty by byl poněkud složitější a obsáhlejší. Příklad 1.2 Metodou regula falsi určete kořen funkce f (x) = cos(x) + x/2 − 1 ležící v intervalu h0,6; 1,5i. Podle věty 1.1, lze snadno ověřit, že v intervalu h0,6; 1,5i leží právě jeden kořen. Jednotlivé aproximace určené metodou regula falsi (viz vztah (1.17) jsou uspořádány do tabulky 2. Postup je patrný z obrázku 5.
Obr. 5: Metoda regula falsi; f (x) = cos(x) + x/2 − 1; x1 = 1,5; x2 = 0,6
Metody řešení nelineárních rovnic i 1 2 3 4 5 6 7 8
xi
16 f (xi ) −0,179 263 0,125 336 0,050 192 0,008 954 0,001 289 0,000 179 0,000 025 3,4.10−6
1,5 0,6 0,970 330 1,086 193 1,105 879 1,108 691 1,109 082 1,109 136
Tab. 2: Metoda regula falsi; f (x) = cos(x) + x/2 − 1; x1 = 1,5; x2 = 0,6
Osmá aproximace kořene ξ je x8 = 1,109 136 s funkční hodnotou f (x8 ) = 3,4.10−6 . Příklad 1.3 Metodou regula falsi určete kořen funkce f (x) = tervalu h0,2; 5i.
1 x2
−
1 2
ležící v in-
Snadno se ověří, že předpoklady věty 1.1 jsou splněny, a tedy, že v intervalu h0,2; 5i leží právě jeden kořen. Kořen ξ můžeme v tomto případě snadno určit algebraicky, tj. √ . ξ = 2 = 1,414 213 562. Posloupnost {xi }∞ i=1 , viz tabulka 3, generovaná metodou regula falsi, konverguje velice pomalu. To je patrné i z obrázku 6. Detailní výřez, vytečkovanou obdélníkovou oblast, zachycuje obrázek 7.
i
xi
i
xi
i
xi
1 0,2
21 3,611 467
101 1,602 202
2 5
22 3,553 522
102 1,595 909
3 4,911 538
23 3,496 895
103 1,589 818
4 4,824 977
24 3,441 561
104 1,583 924
5 4,740 278
25 3,387 494
105 1,578 219
Tab. 3: Metoda regula falsi; f (x) =
1 x2
− 12 ; x1 = 0,2; x2 = 5
Metody řešení nelineárních rovnic
Obr. 6: Metoda regula falsi; f (x) =
17
1 x2
− 12 ; x1 = 0,2; x2 = 5
Obr. 7: Metoda regula falsi; detail obrázku 6
Metody řešení nelineárních rovnic
1.4
18
Metoda sečen
Modifikací předchozí metody je metoda sečen. Zde upouštíme od požadavku, aby funkce f (x) měla různá znaménka v bodech, které používáme pro určení další aproximace. Na rozdíl od předchozích metod, metoda sečen užívá k určení další iterace vždy dvou bezprostředně předcházejících iterací. Říkáme, že metoda sečen je dvoubodová stacionární metoda (viz kap. 1.1 - vztah (1.5)), tj. xi+1 = xi −
xi − xi−1 f (xi ), i = 2, 3, . . . . f (xi ) − f (xi−1 )
(1.21)
Tato metoda znamená oproti metodě regula falsi podstatné a to přede√ zlepšení, . vším v rychlosti konvergence. Řád metody sečen je (1 + 5)/2 = 1,618, viz [6]. Nevýhodou je, že metoda při nevhodné volbě počátečních aproximací nemusí konvergovat. Takový případ je popsán v příkladu 1.5. Příklad 1.4 Metodou sečen určete kořen funkce f (x) = cos(x) + x/2 − 1. Počáteční aproximace zvolte x1 = 1,5; x2 = 0,6. Z příkladu 1.2 už víme, že v intervalu h0,6; 1,5i leží právě jeden kořen. Postup určování aproximací je patrný z obrázku 8. Jednotlivé aproximace jsou zapsány do tabulky 4.
Obr. 8: Metoda sečen; f (x) = cos(x) + x/2 − 1; x1 = 1,5; x2 = 0,6
Metody řešení nelineárních rovnic i 1 2 3 4 5 6 7 8
xi 1,5 0,6 0,970 330 1,217 693 1,100 291 1,108 664 1,109 147 1,109 144
19 f (xi ) −0,179 263 0,125 336 0,050 192 −0,045 342 0,003 482 0,000 190 −9,6.10−7 2,6.10−10
Tab. 4: Metoda sečen; f (x) = cos(x) + x/2 − 1; x1 = 1,5; x2 = 0,6
Získali jsme aproximaci x8 = 1,109 144 s funkční hodnotou f (x8 ) = 2,6.10−10 . Další příklad demonstruje špatnou volbu počátečních aproximací. Příklad 1.5 . Metodou sečen určete kořen funkce f (x) = cos(x) + x/2 − 1 v intervalu h0,3; 1,5i. Počáteční aproximace zvolte x1 = 1,5; x2 = 0,3 (tj. krajní body intervalu). Jelikož funkce f nabývá v krajních bodech intervalu h0,3; 1,5i hodnot s opačnými . . znaménky (f (0,3) = 0,1053 a f (1,5) = −0,1793), leží v tomto intervalu alespoň jeden kořen (viz věta 1.1). Postup při výpočtu dalších aproximací je zřejmý z obrázku 9. Hodnoty těchto aproximací jsou uvedeny v tabulce 5.
i
xi
f (xi )
1
1,5
−0,179 263
2
0,3
0,105 336
3
0,744 147
0,107 740
4 −19,169 034 −9,635 117 5
0,523 940
0,127 825
6
0,266 103
0,097 855
7 8
−0,575 753 −0,449 094 0,115 487
0,051 082
Tab. 5: Metoda sečen; f (x) = cos(x) + x/2 − 1; x1 = 1,5; x2 = 0,3
Metody řešení nelineárních rovnic
20
Obr. 9: Metoda sečen; f (x) = cos(x) + x/2 − 1; x1 = 1,5; x2 = 0,3
Setkáváme se zde se situací, kdy dvě po sobě následující aproximace (v našem případě jsou to x2 a x3 ) mají relativně blízké funkční hodnoty. To má za důsledek, že následující aproximace (x4 ) leží vně intervalu h0,3; 1,5i. Posloupnost {xi }∞ i=1 potom nekonverguje k hledanému kořenu (tj. z intervalu h0,3; 1,5i). Poznámka 1.13 Samozřejmě, může také nastat případ, že se funkční hodnoty dvou po sobě následujících aproximací budou shodovat, tj. f (xi ) = f (xi−1 ), potom další aproximace není dokonce definována, viz vztah (1.21).
1.5
Metoda Newtonova
Metoda sečen je dvoubodovou iterační metodou. Je založena na přímce procházející body [xi , f (xi )] a [xi−1 , f (xi−1 )]. Následující aproximace xi+1 je definována jako nulový bod této přímky. Naproti tomu, Newtonova metoda následující aproximaci xi+1 definuje jako nulový bod přímky, jež je tečnou ke grafu funkce f s dotykem v bodě [xi , f (xi )], viz obrázek 10.
Metody řešení nelineárních rovnic
21
Obr. 10: Newtonova metoda; f (x) = cos(x) + x/2 − 1; x1 = 0,7
Z uvedeného je zřejmé, že Newtonova metoda je jednobodovou iterační metodou. Je dána vztahem xi+1 = xi −
f (xi ) , i = 1, 2, . . . ; f 0 (xi )
(1.22)
tomu odpovídá Newtonova iterační funkce tvaru F (x) = x −
f (x) . f 0 (x)
(1.23)
Poznámka 1.14 Je přirozené, že Newtonově metodě se také často říká metoda tečen. Věta 1.15 Nechť f (x) je funkce spojitá na intervalu ha, bi a nechť na tomto intervalu jsou spojité i její první a druhá derivace, f 0 (x) a f 00 (x). Nechť ξ je kořenem rovnice f (x) = 0 a f 0 (ξ) 6= 0. Pak existuje δ > 0 tak, že posloupnost {xi }∞ i=1 generovaná Newtonovou metodou (1.22) konverguje k bodu ξ pro každou počáteční aproximaci x1 ∈ hξ − δ, ξ + δi ⊆ ha, bi. Důkaz: Ukážeme, že existuje interval hξ − δ, ξ + δi ⊆ ha, bi, na kterém iterační funkce splňuje předpoklady věty 1.6. Protože f 0 (x) je spojitá funkce na intervalu ha, bi a f 0 (ξ) 6= 0, ∃δ1 > 0 : ∀x ∈ hξ − δ1 , ξ + δ1 i : f 0 (x) 6= 0. Funkce F (x) je tedy
Metody řešení nelineárních rovnic
22
definovaná a spojitá na intervalu hξ − δ1 , ξ + δ1 i. Derivací vztahu (1.23) dostaneme F 0 (x) = 1 −
f 02 (x) − f (x)f 00 (x) f (x)f 00 (x) = , pro x ∈ hξ − δ1 , ξ + δ1 i. (1.24) f 02 (x) f 02 (x)
Protože funkce f (x) i její derivace f 0 (x) a f 00 (x) jsou spojité na intervalu ha, bi, jsou funkce F (x) a její derivace F 0 (x) spojité na intervalu hξ − δ1 , ξ + δ1 i. Jelikož f (ξ) = 0, je f (ξ)f 00 (ξ) F 0 (ξ) = = 0. (1.25) f 02 (ξ) Z faktů, že F 0 (ξ) = 0 a F 0 (x) je spojitá funkce na intervalu hξ − δ1 , ξ + δ1 i, plyne, že ∃ δ, 0 < δ < δ1 , takové, že |F 0 (x)| ≤ q < 1 pro každé x ∈ hξ − δ, ξ + δi. Zbývá ukázat, že funkce F splňuje Lipschitzovu podmínku. Ta plyne z věty o střední hodnotě, neboť pro libovolná čísla x, y ∈ hξ − δ, ξ + δi platí |F (x) − F (y)| = |F 0 (α)| |x − y| ≤ q |x − y| , kde α je vhodný bod mezi body x a y. Funkce F (x) splňuje na intervalu hξ − δ, ξ + δi předpoklady věty 1.6, a tudíž posloupnost {xi }∞ i=1 generovaná Newtonovou metodou konverguje k bodu ξ pro každou počáteční aproximaci x1 ∈ hξ − δ, ξ + δi ⊆ ha, bi. 2
Věta 1.16 Newtonova metoda (1.22) je metoda druhého řádu. Důkaz: Důkaz plyne z věty 1.5. Platí totiž • ξ = F (ξ), neboť ξ je pevný bod funkce F ; • F 0 (ξ) = 0, viz vztah (1.25); f 00 (ξ) • F (ξ) = 0 6 0, čož plyne derivováním vztahu (1.24). = f (ξ) 00
2
Příklad 1.6 Newtonovou metodou určete kořen funkce f (x) = cos(x) + x/2 − 1. Počáteční aproximaci zvolte x1 = 0,7. Postup určování aproximací je patrný z obrázku 10. Jednotlivé aproximace jsou zapsány do tabulky 6.
Metody řešení nelineárních rovnic i 1 2 3 4 5 6
23
xi
f (xi ) 0,114 842 −0,177 428 −0,012 191 −0,000 190 −5,2.10−8 −3,8.10−15
0,7 1,496 311 392 1,139 476 135 1,109 625 380 1,109 144 312 1,109 144 182
Tab. 6: Newtonova metoda; f (x) = cos(x) + x/2 − 1; x1 = 0,7
Dospěli jsme k aproximaci x6 f (x6 ) = −3,8.10−15 .
1.6
=
1,109 144 182 s funkční hodnotou
Rozšířená Newtonova metoda
Je-li ξ nulový bod funkce f (x) a existují-li ve všech bodech okolí U (ξ) bodu ξ všechny derivace f (i) (x), které budeme potřebovat, pak Taylorův rozvoj funkce f (x) v bodě x1 ∈ U (ξ) je f (ξ) = 0 = f (x1 ) + (ξ − x1 )f 0 (x1 ) +
(ξ − x1 )2 00 f (x1 ) + . . . . 2!
(1.26)
Jestliže bod x1 je dostatečně blízko kořene ξ, jsou vyšší mocniny (ξ − x1 )i zanedbatelné a rovnici (1.26) lze potom psát ve tvaru např. 0 = f (x1 ) + (x2 − x1 )f 0 (x1 ),
(1.27)
nebo
(x2 − x1 )2 00 f (x1 ). (1.28) 2! Zde místo ξ píšeme x2 , neboť x2 je pouze přiblížením k požadovanému nulovému bodu ξ. Řešením těchto rovnic získáme 0 = f (x1 ) + (x2 − x1 )f 0 (x1 ) +
x2 = x1 − a
f 0 (x1 ) x2 = x1 − 00 ± f (x1 )
p
f (x1 ) f 0 (x1 )
[f 0 (x1 )]2 − 2f (x1 )f 00 (x1 ) , f 00 (x1 )
(1.29)
(1.30)
kde ze dvou kořenů (1.30) rovnice (1.28) je přiblížením ke kořenu ξ rovnice (1.1) právě ten kořen, který má u třetího sčítance ve vztahu (1.30) znaménko souhlasné se znaménkem derivace f 0 (x1 ) (viz poznámka 1.18). Dospěli jsme k iteračním metodám xi+1 = F (xi ), kde iterační funkce jsou F (x) = x −
f (x) f 0 (x)
(1.31)
Metody řešení nelineárních rovnic
24
a
p [f 0 (x)]2 − 2f (x)f 00 (x) f 0 (x) F (x) = x − 00 + sign[f 0 (x)] . (1.32) f (x) f 00 (x) Vidíme, že v prvním případě se jedná o klasickou Newtonovu metodu, tj. metodu tečen - viz kap. 1.5, vztah (1.23). Skutečně, xi+1 je nulovým bodem lineární funkce (přímky), která má v bodě xi stejnou funkční hodnotu a stejnou hodnotu 1. derivace jako funkce f (x), tj. daná přímka je tečnou ke grafu funkce f (x) v bodě xi . Druhý případ je zřejmým rozšířením Newtonovy metody. Aproximace xi+1 je nulovým bodem kvadratické funkce (paraboly), která má v bodě xi stejnou funkční hodnotu a stejné hodnoty 1. a 2. derivace jako funkce f (x). Metodu, jejíž iterační funkcí je funkce (1.32), nazýváme rozšířenou Newtonovou metodou. Poznámka 1.17 Nechť ξ je jednoduchý kořen. Ukážeme, že pokud je reálný bod x dostatečně blízko kořene ξ, nabývá iterační funkce (1.32) reálné hodnoty. Funkce f (x) i f 0 (x) jsou spojité, f (ξ) = 0, f 0 (ξ) 6= 0, proto ∃δ > 0 : ∀x ∈ hξ − δ; ξ + δi : [f 0 (x)]2 − 2f (x)f 00 (x) > 0.
Poznámka 1.18 Bod ξ je pevným bodem iterační funkce F (x), tzn. ξ = F (ξ). Jelikož f (ξ) = 0, ze vztahu (1.32) dostáváme ξ = F (ξ) = ξ −
f 0 (ξ) | f 0 (ξ) | ± . f 00 (ξ) f 00 (ξ)
Je zřejmé, že druhý a třetí člen v předchozím vztahu se musí navzájem odečíst, a tedy ze znamének „±ÿ skutečně volíme to, které je dáno funkcí sign[f 0 (x)]. Poznámka 1.19 Podobně, jako pro Newtonovu metodu požadujeme f 0 (ξ) 6= 0 (viz věta 1.15), požadujeme v případě rozšířené Newtonovy metody f 00 (ξ) 6= 0. Potom pro každou aproximaci xi (dostatečně blízkou kořenu ξ), je aproximace xi+1 = F (xi ), daná iterační funkcí (1.32), definovaná. Zobecnění těchto metod lze získat pomocí Taylorova rozvoje zanedbáním členů (ξ − x1 )k (k) po f (x1 ), k = 1, 2, . . . . Zjištěním první nenulové derivace iterační k! funkce F (x) (viz (1.31) a (1.32)) v bodě ξ lze dospět (na základě věty 1.5) k řádu konvergence dané metody. U klasické Newtonovy metody (1. případ) platí: F 0 (ξ) = 0, F 00 (ξ) 6= 0, a tedy řád klasické Newtonovy metody je 2 (metoda konverguje kvadraticky). U rozšířené metody (2. případ) platí: F 0 (ξ) = 0, F 00 (ξ) = 0, F 000 (ξ) 6= 0, a tedy řád je 3 (metoda konverguje kubicky, tj. rychleji než klasická Newtonova metoda ).
Metody řešení nelineárních rovnic
25
Příklad 1.7 Užitím rozšířené Newtonovy metody určete kořen funkce f (x) = cos(x) + x/2 − 1. Počáteční aproximaci zvolte x1 = 2,5. Proces výpočtu jednotlivých aproximací vidíme na obrázku 11. Jednotlivé aproximace jsou zapsány do tabulky 7.
Obr. 11: Rozšířená Newtonova metoda; f (x) = cos(x) + x/2 − 1; x1 = 2,5
i 1 2 3 4 5
xi 2,5 1,443 507 781 1,122 644 623 1,109 145 115 1,109 144 182
f (xi ) −0,551 144 −0,151 301 −0,005 377 −3,7.10−7 < 10−15
Tab. 7: Rozšířená Newtonova metoda; f (x) = cos(x) + x/2 − 1; x1 = 2,5
Již v páté aproximaci jsme dosáhli hodnoty x5 = 1,109 144 182 s hodnotou funkce f (x5 ) < 10−15 .
Metody řešení nelineárních rovnic
1.7
26
Srovnání metod
Porovnejme jednotlivé metody z hlediska jejich rychlosti konvergence. Při seznamování se s jednotlivými metodami jsme už uvedli jejich řád konvergence. Zjistěme nyní na našem testovacím příkladu, kolik iterací je potřeba provést k dosažení stanovené přesnosti. V tabulce 8 znovu uvedeme posloupnosti iterací {xi }ki=1 pro jednotlivé metody. Přitom, výpočet zastavíme vždy, když relativní chyba xi − xi−1 xi bude menší než 10−6 . řád x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 ··· x21 x22
1 MP 0,785398163 1,570796327 1,178097245 0,981747704 1,079922475 1,129009860 1,104466167 1,116738014 1,110602090 1,107534129 ······ 1,109144509 1,109143760
1 MRF
1,618 MS
2 NM
1,5 0,6 0,970330404 1,086193055 1,105878748 1,108691394 1,109081629 1,109135544 1,109142989 1,109144017
1,5 0,6 0,970330404 1,217693011 1,100290627 1,108664337 1,109146603 1,109144181 1,109144182
0,7 1,496311392 1,139476135 1,109625380 1,109144312 1,109144182
3 RNM 2,5 1,443507781 1,122644623 1,109145115 1,109144182
Tab. 8: Řešení rovnice cos(x) + x/2 − 1 = 0; dosažení rel. chyby < 10−6
Vidíme, že rychlost konvergence posloupností {xi }∞ i=1 určených jednotlivými metodami pro náš testovací příklad odpovídá řádu konvergence metod. Pro získání kořene s relativní chybou menší než 10−6 jsme museli metodou půlení (MP) provést 22 iterací. U metody regula falsi (MRF) nám stačilo 10 iterací. Ještě lepších výsledků jsme dosáhli v případě použití metod s vyšším řádem konvergence. U metody sečen (MS) to bylo 9 iterací, u Newtonovy metody (NM) už jen 6, a u rozšířené Newtonovy metody (RNM) dokonce jen 5 iterací. Musíme si ovšem uvědomit, že v případě Newtonovy metody potřebujeme v každé iteraci vyčíslit hodnotu funkce a hodnotu její první derivace v bodě xi (tj. f (xi ), f 0 (xi )), kdežto v případě rozšířené Newtonovy metody potřebujeme vyčíslit navíc i hodnotu druhé derivace v bodě xi (tj. f (xi ), f 0 (xi ) a f 00 (xi )).
Výpočet kořenů polynomu
2
27
Výpočet kořenů polynomu
V této kapitole se zaměříme na speciální, ale velice častý, typ nelineárních rovnic, na polynomické rovnice. V dalším výkladu budeme uvažovat polynom p stupně n ≥ 2 s reálnými koeficienty tvaru p(x) = an xn + an−1 xn−1 + · · · + a1 x + a0 .
(2.1)
Budeme hledat kořeny polynomu (2.1), tj. hledat řešení rovnice p(x) = 0.
(2.2)
Je vhodné problém určení všech reálných kořenů polynomu p(x) rozdělit na dva dílčí problémy: 1. nalezení největšího kořene; 2. nalezení ostatních kořenů.
2.1
Určení největšího kořene polynomu
Newtonova metoda, xi+1 = xi −
p(xi ) , p0 (xi )
i = 1, 2, . . . ,
(2.3)
umožňuje při vhodné volbě počáteční hodnoty x1 nalézt kořen polynomu (2.1). Následující věta ukazuje, jak pro polynom s reálnými kořeny získat jeho největší kořen. Věta 2.1 Nechť p(x) je polynom stupně n ≥ 2 s reálnými koeficienty tvaru (2.1). Jestliže všechny kořeny ξj , j = 1, 2, . . . , n, polynomu p(x) jsou reálné, ξ1 ≥ ξ2 ≥ · · · ≥ ξn , pak Newtonova metoda dává ryze klesající posloupnost ∞ {xi }i=1 konvergující ke kořenu ξ1 pro každou počáteční hodnotu x1 > ξ1 . Důkaz: Bez újmy na obecnosti můžeme předpokládat, že p(x1 ) > 0. Protože p(x) nemění znaménko pro x > ξ1 , platí p(x) = an xn + an−1 xn−1 + · · · + a1 x + a0 > 0 ∀x > ξ1 , a proto an > 0. Derivace p0 (x) = nan xn−1 + (n − 1)an−1 xn−2 + · · · + 2a2 x + a1 má n − 1 kořenů αj , j = 1, . . . , n − 1, pro něž podle Rolleovy věty platí: ξ1 ≥ α1 ≥ ξ2 ≥ α2 ≥ · · · ≥ αn−1 ≥ ξn .
Výpočet kořenů polynomu
28
Protože p0 (x) je stupně n − 1 ≥ 1, jsou α1 , α2 , . . . , αn−1 všechny jeho kořeny, a platí p0 (x) > 0 pro každé x > α1 , neboť an > 0. Aplikací Rolleovy věty dále získáme p00 (x) > 0 pro x > α1 ,
(2.4)
neboť an > 0, a podobně, s vědomím, že n ≥ 2, p000 (x) ≥ 0 pro x ≥ α1 .
(2.5)
Pro xi > ξ1 plyne, že xi+1 = xi −
p(xi ) < xi , p0 (xi )
protože p(xi ) > 0, p0 (xi ) > 0. Zbývá ukázat, že xi+1 > ξ1 . Ze vztahů (2.4), (2.5), xi > ξ1 ≥ α1 a Taylorovy věty plyne, že 1 0 = p(ξ1 ) = p(xi ) + (ξ1 − xi )p0 (xi ) + (ξ1 − xi )2 p00 (δ) > p(xi ) + (ξ1 − xi )p0 (xi ), 2 kde ξ1 < δ < xi . Ze vztahu (2.3) platí p(xi ) = p0 (xi )(xi − xi+1 ). Dosazením získáme 0 > p0 (xi )(xi − xi+1 + ξ1 − xi ) = p0 (xi )(ξ1 − xi+1 ), a tedy xi+1 > ξ1 , protože p0 (xi ) > 0.
2
Poznámka 2.2 Určit počáteční hodnotu x1 takovou, aby byla větší než největší kořen, lze pro daný polynom poměrně snadno. Každý reálný kořen ξj , j = 1, . . . , n, polynomu (2.1) totiž musí splňovat následující podmínku a1 an−1 a0 . |ξj | ≤ max , 1 + , . . . , 1 + an an an Podmínek, vymezujících hranice kořenů, existuje celá řada. Další jsou uvedeny např. v knize [1]. Newtonova metoda konverguje kvadraticky (řád této iterační metody je 2). Konverguje tedy ke kořenu zpravidla rychleji než některé jednodušší metody (metoda půlení, metoda regula falsi, metoda sečen, viz kapitola 1). Přesto, jestliže ∞ je počáteční hodnota x1 daleko od kořene, tj. x1 ξ1 , pak posloupnost {xi }i=1 získaná Newtonovou metodou konverguje, zvláště v počátku, velmi pomalu. Skutečně, podle (2.3) je xni + · · · 1 xi+1 = xi − n−1 ≈ xi 1 − , i = 1, 2, . . . , n nxi + · · ·
Výpočet kořenů polynomu
29
takže u polynomů vyšších stupňů je rozdíl mezi xi a xi+1 malý. Tato skutečnost vede k úvaze následující dvojkrokové metody xi+1 = xi − 2
p(xi ) , p0 (xi )
i = 1, 2, . . . ,
(2.6)
místo Newtonovy metody (2.3). Je zde nyní nebezpečí „přeskočeníÿ kořene ξ1 . Jsou možné dva případy, jak se bude chovat posloupnost {xi }∞ i=1 definovaná dvojkrokovou metodou (2.6) pro x 1 > ξ1 : 1. x1 ≥ x2 ≥ · · · ≥ xi ≥ xi+1 ≥ · · · ≥ ξ1 a lim xi = ξ1 , i→∞
2. ∃ i0 ∈ N : p(x1 )p(xi ) > 0 pro 1 ≤ i < i0 a p(x1 )p(xi0 ) < 0, tj. x1 > x2 > · · · > xi0 −1 > ξ1 > xi0 . V prvním případě posloupnost {xi }∞ i=1 konverguje monotónně, a rychleji než Newtonova metoda (2.3), ke kořenu ξ1 . Ve druhém případě je i0 -tý člen posloupnosti {xi }∞ i=1 první, který je menší než kořen ξ1 . Lze ho užít jako počátečního bodu y0 := xi0 Newtonovy metody (2.3) yi+1 = yi −
p(yi ) , p0 (yi )
i = 0, 1, . . . ,
také s předpokladem monotónní konvergence y1 ≥ y2 ≥ · · · ≥ ξ1 a lim yi = ξ1 . Platí totiž následující věta. i→∞
Věta 2.3 Nechť p(x) je reálný polynom stupně n ≥ 2 tvaru (2.1). Nechť všechny kořeny ξ1 ≥ ξ2 ≥ · · · ≥ ξn jsou reálné. Nechť α1 je největší kořen p0 (x), a tedy ξ1 ≥ α1 ≥ ξ2 . Pro n = 2 nechť navíc ξ1 > ξ2 . Pak pro každé z > ξ1 jsou čísla z 0 := z −
p(z) , p0 (z)
y := z − 2
p(z) , p0 (z)
y 0 := y −
p(y) p0 (y)
definovaná a platí: α1 < y, ξ1 ≤ y 0 ≤ z 0 . Obrázek 12 znázorňuje geometrickou interpretaci ”dvojkrokové metody”. Důkaz: I zde bez újmy na obecnosti předpokládejme, že an > 0, a tedy p(z) > 0 pro z > ξ1 . Pro takové hodnoty z definujme veličiny: ∆0 := p(z 0 ), ∆1 := p(z 0 ) − p(y) − (z 0 − y)p0 (y).
Výpočet kořenů polynomu
Obr. 12: Znázornění dvojkrokové metody
Obr. 13: Veličiny ∆0 a ∆1 interpretované jako plochy
30
Výpočet kořenů polynomu
31
Veličiny ∆0 a ∆1 mohou být též interpretovány jako příslušné plochy nad a pod grafem funkce p0 (x), viz obrázek 13, neboť Zz0
0
[p0 (t) − p0 (z)] dt = [p(t) − p0 (z) t]zz = p(z 0 ) − p(z) − p0 (z)(z 0 − z) = p(z 0 ) = ∆0
z
Zz0 a
0
[p0 (t) − p0 (y)] dt = [p(t) − p0 (y) t]zy = p(z 0 ) − p(y) − p0 (y)(z 0 − y) = ∆1 .
y
Protože z 0 − y = z − z 0 > 0 a p0 (x) je konvexní pro x ≥ α1 , dostáváme (viz obrázek 13) ∆1 ≤ ∆0 , jestliže y ≥ α1 , (2.7) přičemž ∆1 = ∆0 právě tehdy, když p0 (x) je lineární funkcí, tj. právě tehdy, když p je stupně 2. Nyní rozlišíme 3 případy, y > ξ1 , y = ξ1 , y < ξ1 . Pro y > ξ1 tvrzení věty plyne z Věty 2.1. Pro y = ξ1 ukážeme zaprvé, že ξ1 > α1 > ξ2 , a tedy, že ξ1 je jednoduchý kořen polynomu p(x). Pokud by totiž ξ1 byl násobný kořen, a tedy y = ξ1 = ξ2 = α1 , pak pro n ≥ 3 je ∆1 < ∆0 (plyne z (2.7)), což je ve sporu s tvrzením ∆1 = p(z 0 ) − p(ξ1 ) − (z 0 − ξ1 )p0 (ξ1 ) = p(z 0 ) = ∆0 . Tedy ξ1 musí být jednoduchý kořen. Odtud α1 < ξ1 = y 0 = y < z 0 a tvrzení věty platí i pro tento případ. Zbývá vyšetřit případ y < ξ1 . Uvažujme, že α1 < y (že tomu tak musí být, ukážeme v poslední části důkazu). Protože p(z) > 0 a ξ2 < α1 < y < ξ1 , je p(y) < 0, p0 (y) > 0, a tedy y 0 je definováno. Protože p(y) = (y − y 0 )p0 (y) a ∆1 ≤ ∆0 , dostáváme ∆0 − ∆1 = p(y) + (z 0 − y)p0 (y) = p0 (y)(z 0 − y 0 ) ≥ 0. Proto z 0 ≥ y 0 . Konečně, z Taylorovy věty plyne 1 0 = p(ξ1 ) = p(y) + (ξ1 − y)p0 (y) + (ξ1 − y)2 p00 (δ), y < δ < ξ1 , 2 a protože p00 (x) ≥ 0 pro x ≥ α1 a p(y) = (y − y 0 )p0 (y), platí 0 ≥ p(y) + (ξ1 − y)p0 (y) = p0 (y)(ξ1 − y 0 ). Protože p0 (y) > 0, dostáváme ξ1 ≤ y 0 . K dokončení důkazu stačí ukázat, že y = y(z) > α1 pro každé z > ξ1 .
(2.8)
Výpočet kořenů polynomu
32
Znovu rozlišíme dva případy: ξ1 > α1 > ξ2 a ξ1 = α1 = ξ2 . Jestliže ξ1 > α1 > ξ2 , pak vztah (2.8) platí vždy, když ξ1 < z < ξ1 + (ξ1 − α1 ). To proto, jelikož Věta 2.1 implikuje z > z 0 ≥ ξ1 , a proto y = z 0 − (z − z 0 ) > ξ1 − (ξ1 − α1 ) = α1 . Tudíž můžeme vybrat z0 splňující y(z0 ) > α1 . Předpokládejme, že existuje z1 > ξ1 s podmínkou y(z1 ) ≤ α1 . Z věty o střední hodnotě spojitých funkcí plyne, že existuje z ∈ hz0 , z1 i s podmínkou y = y(z) = α1 . Ze vztahu (2.7) pro z = z platí: ∆1 = p(z 0 ) − p(y) − (z 0 − y)p0 (y) = p(z 0 ) − p(y) ≤ ∆0 = p(z 0 ),
(2.9)
a proto p(y) = p(α1 ) ≥ 0. Na druhou stranu p(α1 ) < 0, protože ξ1 je jednoduchý kořen, což je spor se změnou znaménka p(x). Vztah (2.8) tedy musí platit pro každé z > ξ1 . Jestliže ξ1 = α1 = ξ2 , uvažujme n ≥ 3. Předpokládejme bez újmy na obecnosti, že p(x) = xn + an−1 xn−1 + · · · + a1 x + a0 . Pak 1 + an−1 + · · · + zan0 p(z) z z z z =z− 0 =z− n−1 an−1 1 a1 = z − p (z) n 1 + n z + · · · + n zn−1 n 0
1 1+O . z
Proto 2z y = y(z) = z + 2(z − z) = z − n 0
1 2 1+O =z 1− + O(1). z n
Protože n ≥ 3, hodnota y(z) roste neomezeně s z → ∞ a znovu vidíme, že ∃z0 > ξ1 : y0 = y(z0 ) > α1 . Jestliže vztah (2.8) není splněn pro všechna z > ξ1 , pak ∃z > ξ1 : y = y(z) = α1 . Opět dospějeme ke sporu, neboť ze vztahu (2.9) bychom obdrželi p(y) = p(α1 ) ≥ 0, což není možné. 2 Takto, kombinací dvojkrokové (2.6) a Newtonovy metody (2.3), lze poměrně efektivně určit největší kořen ξ1 polynomu (2.1). Demonstrujme si na následujícím příkladu výhodu dvojkrokové metody (2.6) proti Newtonově metodě (2.3) při určování největšího kořene. Příklad 2.1 Určete největší kořen polynomu 8 Y p(x) = (x − i) = x8 − 36x7 + 546x6 − 4536x5 + i=1
22449x4 − 67284x3 + 118124x2 − 109584x + 40320. Použijte dvojkrokovou metodu a Newtonovu metodu. Porovnejte obě metody z hlediska rychlosti konvergence. Jako počáteční aproximaci zvolte x1 = 20.
Výpočet kořenů polynomu
33
Uvažovaný polynom má kořeny ξ1 = 8, ξ2 = 7, . . . , ξ8 = 1, tzn. x1 = 20 je větší, než největší kořen polynomu, a tedy užitím obou metod s touto počáteční hodnotou musíme skutečně dospět k hodnotě největšího kořene ξ1 = 8. Jednotlivé aproximace určené Newtonovou metodou, tj. užitím vzorce (2.3), stejně jako aproximace určené dvojkrokovou metodou, tj. dané vzorcem (2.6), jsou zapsány do tabulky 9. xi x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12
Newtonova metoda dvojkroková metoda 20,000 000 20,000 000 18,105 567 16,211 133 16,454 192 13,398 883 15,016 438 11,329 903 13,766 710 9,834 383 12,682 811 8,794 966 11,745 573 8,148 323 10,938 548 7,929 357 10,247 782 8,016 696 9,661 673 8,000 686 9,170 955 8,000 001 8,768 867 8,000 000
Tab. 9: Výpočet největšího kořene Newtonovou metodou a dvojkrokovou metodou
Z tab. 9 je patrné, že posloupnost {xi }∞ i=1 daná dvojkrokovou metodou klesá ke kořenu ξ1 = 8 rychleji než Newtonova metoda, a to až do osmé iterace. U dvojkrokové metody je totiž hodnota x8 oním ”přeskokem” kořene ξ1 = 8. Dál výpočet pokračuje Newtonovou metodou (y0 = x8 ). Můžeme si rovněž povšimnout, že pro následující aproximaci platí ξ1 < x9 < x7 (viz věta 2.3).
2.2
Určení ostatních kořenů polynomu
V minulém odstavci jsme ukázali, jak určit největší kořen ξ1 polynomu (2.1). Zbývá ještě určit ostatní kořeny ξ2 , ξ3 , . . . , ξn polynomu (2.1). Nabízí se dělení kořenovým činitelem x − ξ1 . Uvažujme tedy polynom p1 (x) :=
p(x) . x − ξ1
(2.10)
Největší kořen polynomu p1 (x) je ξ2 a může být určen výše popsaným postupem. Zde ξ1 , nebo ještě lépe y0 = xi0 , může posloužit jako počáteční aproximace. Tímto způsobem lze nakonec získat všechny kořeny. Postupně bychom vyčíslili další polynomy pj−1 (x) pj (x) := , j = 2, 3, . . . , n − 1, (2.11) x − ξj a pro každý potom určili jeho největší kořen ξj+1 .
Výpočet kořenů polynomu
34
Dělení kořenovým činitelem ale není bez rizika, protože zaokrouhlení kořene ξ1 (obecně nejsme schopni vyčíslit ξ1 zcela přesně) znemožní přesně určit polynom p1 (x). Polynom získaný pomocí vztahu (2.10) pro zaokrouhlený kořen ξ1 má kořeny odlišné od ξ2 , ξ3 , . . . , ξn , viz poznámka 2.5. Popsaná chyba se kumuluje a negativně ovlivňuje určování dalších kořenů, viz příklad 2.3. Poznámka 2.4 Ve zbývající části kapitoly budeme užívat označení ξ j , j = 1, 2, . . . , n, pro přibližné hodnoty kořenů ξj , j = 1, 2, . . . , n. Poznámka 2.5 Prozkoumejme podrobněji podíl polynomu p(x) a dvojčlenu x − x0 , kde x0 je libovolné reálné číslo. Platí: p(x) = (x − x0 )qx0 (x) + R,
(2.12)
kde qx0 (x) = bn−1 xn−1 + bn−2 xn−2 + · · · + b1 x + b0 je polynom stupně (n − 1), R = R(x0 ) je zbytek po dělení polynomu p(x) dvojčlenem x − x0 . Číselně je tento zbytek roven hodnotě polynomu p(x) v bodě x0 , tj. R = p(x0 ), což plyne přímo ze vztahu (2.12). Koeficienty bj , j = 0, 1, . . . , n − 1 lze počítat rekurentně pomocí algoritmu známého pod názvem Hornerovo schéma: bn−1 = an br = br+1 x0 + ar+1 ,
r = n − 2, n − 3, . . . , 2, 1, 0.
(2.13)
Pro zbytek R platí R = b−1 = b0 x0 + a0 . Při užití metody dělení kořenovým činitelem potom určujeme druhý největší kořen ξ2 polynomu (2.1) jakožto největší kořen polynomu qξ1 (x) 6= p1 (x). Tyto souvislosti jsou demonstrovány v příkladu 2.2. Při určování dalších kořenů jsou i tyto ovlivněny chybami v určení předchozích kořenů. Příklad 2.2 Uvažujme polynom 4 Y p(x) = (x − i) = x4 − 10x3 + 35x2 − 50x + 24.
(2.14)
i=1
Polynom (2.14) má přesné kořeny ξ1 = 4, ξ2 = 3, ξ3 = 2, ξ4 = 1. Vydělením polynomu (2.14) kořenovým činitelem x − ξ1 = x − 4 dostaneme 3
Y p(x) p1 (x) = = (x − i) = x3 − 6x2 + 11x − 6. x − 4 i=1
(2.15)
Předpokládejme, že největší kořen aproximujeme hodnotou ξ 1 = 4,001. Polynom (2.14) v tomto bodě nabývá hodnoty R = p(ξ 1 ) = 0,006 011 006 001. Pomocí Hornerova schematu (2.13) určíme polynom qξ1 , tj. qξ1 (x) = x3 − 5,999x2 + 10,998 001x − 5,996 997 999.
(2.16)
Největší kořen tohoto polynomu (s přesností 10−9 ) je ξ 2 = 2,996 989 940. Tato hodnota se významně liší od hodnoty ξ2 = 3 druhého největšího kořene polynomu (2.14).
Výpočet kořenů polynomu
35
Problém popsaný výše vyřešil H. Maehly v roce 1954. Vyjádříme-li derivaci polynomu p1 (x) ze vztahu (2.10), tj. 0
p (x) p(x) − , p1 (x) = x − ξ1 (x − ξ1 )2 0
lze psát Newtonovu metodu (2.3) pro polynom p1 (x) ve tvaru xi+1 = xi −
p1 (xi ) p(xi ) , = xi − 0 0 i) p1 (xi ) p (xi ) − (xp(x i −ξ1 )
i = 1, 2, . . . ,
(2.17)
kde není třeba mít polynom p1 (x) vyčíslený. Posloupnost {xi }∞ i=1 daná metodou (2.17) konverguje k přesnému kořenu ξ2 polynomu p(x). Obecně, pro určení kořene ξj+1 zjišťujeme polynomy pj (x) :=
p(x) , (x − ξ1 ) . . . (x − ξj )
j = 1, 2, . . . , n − 1,
a jejich derivace j
0
X p (x) p(x) 1 pj (x) = − . (x − ξ1 ) . . . (x − ξj ) (x − ξ1 ) . . . (x − ξj ) k=1 (x − ξk ) 0
Maehlyho verze Newtonovy metody pro určení kořene ξj+1 , j = 1, 2, . . . , n − 1 je proto následující: xi+1 = xi −
p(xi ) , j P 0 p(xi ) p (xi ) − (xi −ξk )
i = 1, 2, . . . .
(2.18)
k=1
Výhoda vzorce (2.18) spočívá ve skutečnosti, že posloupnost {xi }∞ i=1 konverguje kvadraticky ke kořenu ξj+1 , dokonce i v případě, že čísla ξ1 , ξ2 , . . . , ξj nejsou kořeny polynomu p(x). Tedy výpočet ξj+1 není citlivý na chyby získané ve výpočtu předešlých kořenů. Ukažme na následujícím příkladu výpočet všech kořenů daného polynomu. Příklad 2.3 Určete všechny kořeny polynomu 15 Y p(x) = (x − i) = x15 − 120x14 + . . . , i=1
a to metodou dělení kořenovým činitelem a Maehlyho metodou. V obou případech zvolte jako počáteční aproximaci hodnotu x1 = 20. Nejprve je Newtonovou metodou určen největší kořen polynomu p(x), tj. hodnota ξ 1 . Kořeny ξ j , j = 2, 3, . . . , 15, jsou získány Newtonovou metodou vždy jakožto největší kořen polynomů qξj , j = 1, 2, . . . , 14 (viz poznámka 2.5), resp.
Výpočet kořenů polynomu
36
Maehlyho metodou (2.18). Vždy byl každý kořen zpřesňován tak dlouho, až se dvě po sobě jdoucí aproximace lišily o méně než 10−3 . V tabulce 10 je vidět nepřesnost v kořenech nalezených metodou dělení kořenovým činitelem. ξj ξ1 ξ2 ξ3 ξ4 ξ5 ξ6 ξ7 ξ8 ξ9 ξ 10 ξ 11 ξ 12 ξ 13 ξ 14 ξ 15
dělení koř. činitelem 15, 000 000 036 14, 000 001 944 12, 999 972 698 12, 000 165 123 10, 999 395 686 10, 001 500 466 8, 997 307 628 8, 003 565 277 6, 996 451 559 6, 002 656 420 4, 998 529 021 4, 000 591 721 2, 999 836 728 2, 000 028 914 0, 999 996 781
Maehlyho metoda 15, 000 000 036 14, 000 002 419 13, 000 001 168 12, 000 000 845 10, 999 998 855 10, 000 000 367 8, 999 998 958 8, 000 000 071 6, 999 999 905 6, 000 000 014 5, 000 000 007 4, 000 000 012 3, 000 000 006 2, 000 000 976 1, 000 000 000
Tab. 10: Výpočet kořenů pomocí dělení kořenovým činitelem a Maehlyho metodou
Z uvedeného vyplývá, že spojením dvojkrokové a Maehlyho metody se nabízí velice efektivní způsob výpočtu všech reálných kořenů polynomu s reálnými koeficienty. Algoritmus této kombinované metody lze najít např. v [1].
Citlivost kořenů polynomu
3
37
Citlivost kořenů polynomu
Existuje řada algoritmů pro určení všech kořenů libovolného polynomu. Přesto se může řešitel při výpočtu kořenů polynomu dostat do jistých nesnází. Pokud je například daný polynom určen pouze přibližně, tj. koeficienty polynomu jsou určeny s nějakou nepřesností, byly zaokrouhleny, apod., tato nepřesnost se může mnohonásobně projevit v nepřesnosti hodnot kořenů. Právě citlivost kořenů polynomu na změnu jeho koeficientů je předmětem studia této kapitoly. Podrobněji se tomuto tématu věnují práce [1, 8, 10, 14, 15, 19].
3.1
Špatně podmíněné polynomy
Uvažujme polynom p stupně n ≥ 2 s reálnými koeficienty tvaru p(x) = an xn + an−1 xn−1 + · · · + a1 x + a0 .
(3.1)
Označme kořeny polynomu (3.1) ξj , pro j = 1, 2, . . . , n. Hledáním kořenů polynomu se věnuje předchozí kapitola. Platí tedy p(ξj ) = 0, pro j = 1, 2, . . . , n. Budeme se zabývat tím, jak se změní kořeny polynomu, pokud se polynom (3.1) mírně změní. Budeme určovat kořeny polynomu p(x) + δp(x) ≡ p(x) + δar xr , r = 0, 1, 2, . . . , n ,
(3.2)
tj. uvažujeme pouze změnu koeficientu u mocniny xr , a to z ar na ar + δar . Označíme-li δξj změnu kořene ξj při přechodu od polynomu (3.1) k polynomu (3.2), platí p(ξj + δξj ) + δar (ξj + δξj )r = 0. (3.3) Odtud můžeme vyjádřit p(ξj + δξj ) = −δar (ξj + δξj )r ∼ = −δar ξjr ,
(3.4)
kde přibližný vztah je platný za podmínky δξj ξj . Taylorův rozvoj funkce p(x) v bodě x = ξj je 00
p (ξj ) (δξj )2 + · · · . p(ξj + δξj ) = p(ξj ) + p (ξj )δξj + 2! 0
(3.5)
Platí p(ξj ) = 0. 0 Pokud ξj je jednoduchý kořen, platí také p (ξj ) 6= 0. Jestliže |δξj | 1, můžeme ostatní členy Taylorova rozvoje (3.5) zanedbat. Nyní, odtud porovnáním se vztahem (3.4) získáme vyjádření změny δξj jednoduchého kořene ξj δξj = −
ξjr δar . p0 (ξj )
(3.6)
Citlivost kořenů polynomu
38 0
00
Podobně, pokud ξj je dvojnásobný kořen, platí p (ξj ) = 0 a p (ξj ) 6= 0. I zde uvažujeme pouze první nenulový člen Taylorova rozvoje (3.5), ostatní členy Taylorova rozvoje zanedbáme. Platí tedy (δξj )2 = −
2ξjr δar . p00 (ξj )
(3.7)
Vztahy (3.6) a (3.7) můžeme zobecnit pro k−násobný (k ∈ N) kořen ξj . V tom 0 00 případě je p (ξj ) = p (ξj ) = · · · = p(k−1) (ξj ) = 0 a p(k) (ξj ) 6= 0, a tudíž (δξj )k = −
k!ξjr δar , j = 1, 2, . . . , n a r = 0, 1, . . . , n . p(k) (ξj )
(3.8)
Pokud u nějakého polynomu malá změna koeficientu (změna na vstupu) způsobí mnohem větší změnu jeho kořenů (změna na výstupu), nazýváme takový polynom špatně podmíněným polynomem. Tuto podmíněnost lze definovat jako podíl relativní změny na výstupu a relativní změny na vstupu. K tomu účelu definujme koeficient podmíněnosti kjr jakožto podíl relativní změny kořene ξj (j = 1, 2, . . . , n) a relativní změny koeficientu ar (r = 0, 1, . . . , n), tzn. kjr =
δξj ξj δar ar
pro j = 1, 2, . . . , n a r = 0, 1, . . . , n.
(3.9)
Koeficient podmíněnosti kjr je příslušný j-tému kořenu a r-tému koeficientu. Vezmeme-li v úvahu vztah (3.6), můžeme psát kjr =
ξjr−1 ar p0 (ξj )
pro j = 1, 2, . . . , n a r = 0, 1, . . . , n.
(3.10)
Vztah (3.10) vyjadřuje koeficient podmíněnosti pro jednoduchý kořen ξj . V případě k−násobného kořene ξj využijeme vztah (3.8). Dostáváme 1
kjr =
r
k! k ξjk
−1 1 k 1 δar
[p(k) (ξj )] k
−1
ar
pro j = 1, 2, . . . , n a r = 0, 1, . . . , n.
(3.11)
Číslo podmíněnosti Cp polynomu definujeme jako maximum z absolutních hodnot všech koeficientů podmíněnosti, tj. Cp = max |kjr |. j,r
(3.12)
Pomocí čísla podmíněnosti můžeme definovat špatně podmíněný polynom jako polynom, jehož číslo podmíněnosti je mnohem větší než jedna, tj. Cp 1.
Citlivost kořenů polynomu
39
Poznámka 3.1 Prozkoumejme podrobněji vztah (3.11) udávající hodnotu koeficientu podmíněnosti kjr pro k−násobný kořen ξj . Připomeňme, že koeficient podmíněnosti kjr je definovaný vztahem (3.9) pro (nekonečně) malou změnu δar koeficientu ar , r = 0, 1, . . . , n , polynomu (3.1). Zatímco pro jednoduchý kořen ξj koeficient podmíněnosti kjr nezávisí na velikosti změny δar koeficientu ar (viz vztah (3.10)), u vícenásobného kořene tomu tak není. Ze vztahu (3.11) plyne, že koeficient podmíněnosti kjr k−násobného 1
kořene ξj je přímo úměrný veličině δark nebo rovna dvěma, tj. k ≥ 2, potom
−1
. Je-li tedy násobnost k kořene ξj větší
kjr −→ ±∞ pro δar −→ 0, a tedy také Cp −→ ∞ pro δar −→ 0, viz také příklad 3.4. Z toho vyplývá, že každý polynom mající kořen násobnosti vyšší než jedna je (podle definničního vztahu (3.9)) špatně podmíněným polynomem. U takových polynomů nemá tudíž smysl určovat hodnotu čísla podmíněnosti Cp . Proto se otázkou, zda daný polynom je či není špatně podmíněný budeme zabývat pouze u polynomů, jehož všechny kořeny jsou jednoduché. Poznámka 3.2 Vztahem (3.9) jsme definovali koeficient podmíněnosti kjr . Vyδξ r jadřuje podíl relativní změny ξjj kořene ξj (j = 1, 2, . . . , n) a relativní změny δa ar koeficientu ar (r = 0, 1, . . . , n). Vyjádřeme ještě podíl absolutní změny δξj kořene ξj a absolutní změny δar koeficientu ar . Tento podíl budeme nazývat absolutní koeficient podmíněnosti a označovat Cjr . Tedy Cjr =
δξj δar
pro j = 1, 2, . . . , n a r = 0, 1, . . . , n.
(3.13)
Stejně jako koeficient podmíněnosti kjr , tak i absolutní koeficient podmíněnosti Cjr má smysl určovat pouze v případě jednoduchého kořene (viz poznámka 3.1). Pro něj, ze vztahu (3.6), platí ξjr δξj Cjr = =− 0 δar p (ξj )
3.2
pro j = 1, 2, . . . , n a r = 0, 1, . . . , n.
(3.14)
Příklady polynomů druhého stupně
Zkoumejme nejprve u polynomů druhého stupně, jak se změna některého z koeficientů polynomu projeví na změně jeho kořenů. Příklad 3.1 Uvažujme jednoduchý polynom p(x) = x2 − 3x + 2. Jeho kořeny jsou ξ1 = 1 a ξ2 = 2, viz obr. 14.
(3.15)
Citlivost kořenů polynomu
40
Prozkoumejme, jak bude vypadat polynom (3.15) a jeho kořeny, pokud změníme postupně všechny jeho koeficienty o hodnotu +0,1 resp. −0,1. Na obrázku 15 jsou grafy polynomů 2 − 2 p+ 0 (x) = x − 3x + 2,1 a p0 (x) = x − 3x + 1,9.
Obrázek 16 znázorňuje polynomy 2 2 − p+ 1 (x) = x − 2,9x + 2 a p1 (x) = x − 3,1x + 2
a obrázek 17 zobrazuje polynomy 2 − 2 p+ 2 (x) = 1,1x − 3x + 2 a p2 (x) = 0,9x − 3x + 2.
Obr. 14: p(x) = x2 − 3x + 2
Obr. 15: δa0 = ±0,1
Obr. 16: δa1 = ±0,1
Obr. 17: δa2 = ±0,1
Vypočítejme kořeny takových polynomů a odchylky těchto kořenů od kořenů původního polynomu (3.15). Jejich hodnoty (zaokrouhlené) u všech těchto polynomů uvádí tabulka 11.
Citlivost kořenů polynomu p(x) x − 3x + 2 x2 − 3x + 2, 1 x2 − 3x + 1, 9 x2 − 2, 9x + 2 x2 − 3, 1x + 2 1, 1x2 − 3x + 2 0, 9x2 − 3x + 2 2
ξ1 1 1, 1127 0, 9084 1, 1298 0, 9156 1, 1604 0, 9213
41 ξ2 δξ1 δξ2 2 1, 8873 0, 1127 −0, 1127 2, 0916 −0, 0916 0, 0916 1, 7702 0, 1298 −0, 2298 2, 1844 −0, 0844 0, 1844 1, 5669 0, 1604 −0, 4331 2, 4120 −0, 0787 0, 4120
Tab. 11: Hodnoty kořenů a změn kořenů pro δar = ±0,1
Podobně můžeme uvažovat i změnu polynomu (3.15) např. desetkrát menší, tj. změnu koeficientu ar (r = 0, 1, 2) o hodnotu +0,01; resp. −0,01; viz tabulka 12.
p(x) x − 3x + 2 2 x − 3x + 2, 01 x2 − 3x + 1, 99 x2 − 2, 99x + 2 x2 − 3, 01x + 2 1, 01x2 − 3x + 2 0, 99x2 − 3x + 2 2
ξ1 1 1, 0101 0, 9901 1, 0102 0, 9902 1, 0103 0, 9903
ξ2 δξ1 δξ2 2 1, 9899 0, 0101 −0, 0101 2, 0099 −0, 0089 0, 0099 1, 9798 0, 0102 −0, 0202 2, 0198 −0, 0098 0, 0198 1, 9600 0, 0103 −0, 0400 2, 0400 −0, 0097 0, 0400
Tab. 12: Hodnoty kořenů a změn kořenů pro δar = ±0,01
Aplikací vzorce (3.6) pro polynom (3.15) získáme ξjr δξj = − δar , 2ξj − 3 Konkrétně δξ1 = −
pro j = 1, 2, r = 0, 1, 2.
1r δar = δar , 2·1−3
pro r = 0, 1, 2 ; −δa0 r 2 −2δa1 . δξ2 = − δar = −2r δar = 2·2−3 −4δa2
(3.16)
(3.17)
(3.18)
Pro změnu δa0 = ±0,1 tak dostáváme změny kořenů δξ1 = ±0,1 a δξ2 = ∓0,1. Podobně pro δa1 = ±0,1 je δξ1 = ±0,1 a δξ2 = ∓0,2; resp. pro δa2 = ±0,1 je δξ1 = ±0,1 a δξ2 = ∓0,4. Změně δar = ±0,01 odpovídá změna kořenů úměrně menší. Není bez zajímavosti srovnání těchto přibližných výsledků (jejich první aproximace) se skutečnými, viz srovnání výsledků (3.17) a (3.18) pro δar = ±0,1; resp. δar = ±0,01; s údaji v tabulce 11, resp. s údaji v tabulce 12.
Citlivost kořenů polynomu
42
Ze vztahu (3.10) můžeme určit koeficienty podmíněnosti kjr pro j = 1, 2, r = 0, 1, 2 (viz tabulka 13). j/r 1 2
0 1 2 −2 3 −1 1 −3 2
Tab. 13: Koeficienty podmíněnosti kjr pro polynom p(x) = x2 − 3x + 2
Číslo podmíněnosti je potom Cp = 3, tzn. polynom (3.15) je dobře podmíněným polynomem. Příklad 3.2 Mějme polynom p(x) = x2 − 201x + 10100,
(3.19)
viz obrázek 18. Jeho kořeny, ξ1 = 100 a ξ2 = 101, se liší opět o jedničku (jako v příkladu 1), jen jsou mnohem větší.
Obr. 18: p(x) = x2 − 201x + 10100
Obr. 19: δa0 = ±10−1
Obr. 20: δa1 = ±10−3
Obr. 21: δa2 = ±10−5
Citlivost kořenů polynomu
43
Prozkoumejme, jak bude vypadat polynom (3.19) a jeho kořeny, pokud postupně nepatrně změníme všechny jeho koeficienty. Tyto změněné polynomy jsou graficky znázorněné na obrázcích 19, 20 a 21. V tabulce 14 jsou podle vztahu (3.10) vypočteny koeficienty podmíněnosti kjr . Z tabulky můžeme zjistit hodnotu čísla podmíněnosti polynomu (3.19), Cp = 201, a tudíž polynom (3.19) považujeme za špatně podmíněný polynom. j/r 1 2
0 1 2 −101 201 −100 100 −201 101
Tab. 14: Koeficienty podmíněnosti kjr pro polynom p(x) = x2 − 201x + 10100
Příklad 3.3 Zabývejme se polynomem p(x) = x2 − 3,01x + 2,265.
(3.20)
Jeho kořeny, ξ1 = 1,5 a ξ2 = 1,51, jsou narozdíl od příkladu 1 vzájemně mnohem bližší. Koeficienty podmíněnosti kjr jsou uvedeny v tabulce 15 (vypočteny podle vztahu (3.10)). j/r 1 2
0 1 2 −151 301 −150 150 −301 151
Tab. 15: Koeficienty podmíněnosti kjr pro polynom p(x) = x2 − 3, 01x + 2, 265
Z tabulky můžeme odečíst hodnotu čísla podmíněnosti polynomu (3.20), Cp = 301, a tedy polynom (3.20) je špatně podmíněným polynomem. Příklad 3.4 Uvažujme nyní polynom p(x) = (x − 1,5)2 = x2 − 3x + 2,25.
(3.21)
Polynom (3.21) má dvojnásobný kořen ξ1 = ξ2 = 1,5 ; viz obr. 22. Platí: p0 (x) = 2x − 3 ⇒ p0 (ξj ) = 0, pro j = 1, 2 ; p00 (x) = 2 ⇒ p00 (ξj ) = 2, pro j = 1, 2 . Prozkoumejme, jak bude vypadat polynom (3.21) a jeho kořeny, pokud změníme postupně všechny jeho koeficienty o hodnotu +0,01; resp. −0,01. Na obrázku 23 jsou grafy polynomů 2 − 2 p+ 0 (x) = x − 3x + 2,26 a p0 (x) = x − 3x + 2,24.
Citlivost kořenů polynomu
44
Obrázek 24 znázorňuje polynomy 2 − 2 p+ 1 (x) = x − 2,99x + 2,25 a p1 (x) = x − 3,01x + 2,25
a obrázek 25 zobrazuje polynomy 2 − 2 p+ 2 (x) = 1,01x − 3x + 2,25 a p2 (x) = 0,99x − 3x + 2,25.
Obr. 22: p(x) = x2 − 3x + 2,25
Obr. 23: δa0 = ±0,01
Obr. 24: δa1 = ±0,01
Obr. 25: δa2 = ±0,01
Vypočteme kořeny ξ1 a ξ2 takových polynomů a odchylky δξ1 a δξ2 těchto kořenů od původních kořenů. Jejich hodnoty u všech těchto polynomů uvádí tabulka 16.
Citlivost kořenů polynomu p(x) x − 3x + 2,25 x2 − 3x + 2,26 x2 − 3x + 2,24 x2 − 2,99x + 2,25 x2 − 3,01x + 2,25 1,01x2 − 3x + 2,25 0,99x2 − 3x + 2,25 2
45
ξ1
ξ2
δξ1
δξ2
1,5 1,500 − 0,100i 1,400 1,495 − 0,122i 1,382 1,485 − 0,148i 1,364
1,5 1,500 + 0,100i 1,600 1,495 + 0,122i 1,628 1,485 + 0,148i 1,667
−0,100i −0,100 −0,005 − 0,122i −0,118 −0,015 − 0,148i −0,136
0,100i 0,100 −0,005 + 0,122i 0,128 −0,015 + 0,148i 0,167
Tab. 16: Hodnoty kořenů a změn kořenů pro δar = ±0,01
Skutečné (byť zaokrouhlené) hodnoty δξ1 a δξ2 (z tabulky 16) můžeme porovnat s přibližnými výsledky, které získáme pomocí vztahu (3.7), tzn. s 2ξjr δξ1 = − − 00 δar , p (ξ1 ) s δξ2 =
−
2ξjr δar . p00 (ξ2 )
Tyto přibližné změny δξ1 a δξ2 kořenů polynomu (3.21) pro δar = ±0,01 (r = 0, 1, 2) jsou uvedeny v tabulce 17. p(x) x − 3x + 2,26 x2 − 3x + 2,24 2 x − 2,99x + 2,25 x2 − 3,01x + 2,25 1,01x2 − 3x + 2,25 0,99x2 − 3x + 2,25 2
δξ1 −0,100i −0,100 −0,122i −0,122 −0,150i −0,150
δξ2 0,100i 0,100 0,122i 0,122 0,150i 0,150
Tab. 17: Přibližné hodnoty změn kořenů pro δar = ±0,01
Je zřejmé, že dvojnásobný kořen se změnou koeficientu ar , r = 0, 1, 2 , polynomu (3.21) rozštěpí buď na dva různé reálné kořeny (δar < 0) nebo na dva komplexně sdružené kořeny (δar > 0). Ze vztahu (3.11) můžeme vyjádřit koeficient podmíněnosti r
1
kjr =
2 2 1,5 2 −1 1
22
1 −1 2
δar
r
1,5 2 −1 ar ar = √ δar
pro j = 1, 2 a r = 0, 1, 2.
Hodnoty koeficientů podmíněnosti kjr pro δar = ±0,01 uvádí tabulka 18.
(3.22)
Citlivost kořenů polynomu
δa0 δa0 δa1 δa1 δa2 δa2
46
δar = 0,01 =−0,01 = 0,01 =−0,01 = 0,01 =−0,01
k1r −15i −15 √ −20 1,5i √ −20 1,5 −10i −10
k2r 15i 15 √ 20 1,5i √ 20 1,5 10i 10
Tab. 18: Koeficienty podmíněnosti kjr pro δar = ±0,01 pro polynom (3.21)
Ze vztahů (3.22) i (3.11) je ale zřejmé, že čím menší změnu δar uvažujeme, tím vyšší jsou hodnoty koeficientů podmíněnosti kjr (j = 1, 2 , r = 0, 1, 2). Pro δar −→ 0 tak každý z koeficientů podmíněnosti kjr a tedy také číslo podmíněnosti Cp roste nade všechny meze. To je ovšem zcela v souladu s tvrzením v poznámce 3.1.
3.3
Polynom vyššího stupně
V následujícím příkladu budeme studovat podmíněnost polynomu osmého stupně. Příklad 3.5 Mějme polynom 8 Y p(x) = (x − i) = x8 − 36x7 + 546x6 − 4536x5 + 22449x4 + i=1
−67284x3 + 118124x2 − 109584x + 40320. (3.23) Je to polynom osmého stupně s kořeny ξ1 = 1, ξ2 = 2, . . . , ξ8 = 8, viz obr. 26. Pro vyjádření změn δξj (viz vztah (3.6)) kořenů ξj , j = 1, 2, . . . , 8, polynomu (3.23) si nejdříve vyjádříme hodnoty p0 (ξj ), j = 1, 2, . . . , 8 . Obecně pro každý polynom s kořeny ξ1 , ξ2 , . . . , ξn platí ! n X Y p0 (x) = (x − ξi ) , ∀x ∈ R. j=1
i6=j
Proto p0 (ξj ) =
Y (ξj − ξi ),
pro j = 1, 2, . . . , n ,
i6=j
a konkrétně pro polynom (3.23) dostáváme p0 (ξ) = (−1)ξ (ξ − 1)!(8 − ξ)!
pro ξ = 1, 2, . . . , 8 .
Proto vztah (3.6) dává pro kořen ξj , j = 1, 2, . . . , 8, polynomu (3.23) změnu δξj = (−1)ξj −1
ξjr δar . (ξj − 1)!(8 − ξj )!
(3.24)
Citlivost kořenů polynomu
47
Obr. 26: Polynomy p, p1 , p2
Odtud, anebo přímo ze vztahu (3.14), můžeme vyjádřit hodnoty absolutních koeficientů podmíněnosti Cjr =
ξjr δξj = (−1)ξj −1 δar (ξj − 1)!(8 − ξj )!
Tyto hodnoty jsou obsaženy v tabulce 19. (Pouze z důvodu úspory místa jsou v tabulce zobrazeny absolutní koeficienty podmíněnosti Cjr jen pro r = 0, 1, 3, 5, 7, 8.)
j/r 1 2 3 4 5 6 7 8
0 0,0002 -0,0014 0,0042 -0,0069 0,0069 -0,0042 0,0014 -0,0002
1 0,0002 -0,0028 0,0125 -0,0278 0,0347 -0,0250 0,0097 -0,0016
3 5 7 8 0,0002 0,0002 0,0002 0,0002 -0,0111 -0,0444 -0,1778 -0,3556 0,1125 1,0125 9,1125 27,3375 -0,4444 -7,1111 -113,7778 -455,1111 0,8681 21,7014 542,5347 2712,6736 -0,9000 -32,4000 -1166,4000 -6998,4000 0,4764 23,3431 1143,8097 8006,6681 -0,1016 -6,5016 -416,1016 -3328,8127
Tab. 19: Absolutní koeficienty podmíněnosti Cjr polynomu (3.23)
Citlivost kořenů polynomu
48
Pokud vyčíslíme koeficienty podmíněnosti kjr podle vztahu (3.10), viz tabulka 20, snadno stanovíme číslo podmíněnosti Cp polynomu (3.23), tj. Cp = max |kjr | = |k65 | = 24494, 4 . j,r
j/r 1 2 3 4 5 6 7 8
0 1 3 5 7 8 -8,000 21,743 13,350 0,900 0,007 0,000 28,000 -152,200 -373,800 -100,800 -3,200 0,178 -56,000 456,600 2523,150 1530,900 109,350 -9,113 70,000 -761,000 -7476,000 -8064,000 -1024,000 113,778 -56,000 761,000 11681,250 19687,500 3906,250 -542,535 28,000 -456,600 -10092,600 -24494,400 -6998,400 1166,400 -8,000 152,200 4579,050 15126,300 5882,450 -1143,810 1,000 -21,743 -854,400 -3686,400 -1872,457 416,102 Tab. 20: Koeficienty podmíněnosti kjr polynomu (3.23)
Polynom (3.23) je tedy špatně podmíněným polynomem. Uvažujme nyní dva další polynomy p1 (x) = p(x) + 2.10−4 x7 , p2 (x) = p(x) + 4.10−4 x7 , viz obrázek 26. Vzhledem k polynomu (3.23) se jedná pouze o relativně malou změnu v koefici. . 7 7 entu a7 (a7 = −36, a tedy | δa | = 5,6.10−6 , resp. | δa | = 1,1.10−5 ). a7 a7 Budeme se zajímat o to, jak se tato změna projeví na změně kořenů ξj , j = 1, 2, . . . , 8. j = ξj (p) 1 2 3 4 5 6 7 8
ξj (p1 ) 1,0000 2,0000 3,0018 3,9782 5,1284 5,7714 7,2194 7,9005
ξ j (p1 ) 1,0000 2,0000 3,0018 3,9772 5,1085 5,7667 7,2288 7,9168
ξj (p2 ) 1,0000 1,9999 3,0037 3,9581 5,4141 − 0,1151i 5,4141 + 0,1151i 7,5124 7,6973
ξ j (p2 ) 1,0000 1,9999 3,0036 3,9545 5,2170 5,5334 7,4575 7,8336
Tab. 21: Skutečné hodnoty ξj a aproximované hodnoty ξ j kořenů polynomů p1 a p2
Citlivost kořenů polynomu
49
Tabulka 21 obsahuje aproximované hodnoty ξ j skutečných kořenů ξj polynomů p1 a p2 . Pro tyto hodnoty, viz vztah (3.13), platí ξ j (pi ) = ξj (p) + Cjr δar (pi ),
i = 1, 2, j = 1, 2, . . . , 8, r = 0, 1, 2, . . . , 8 .
Zde ξj (p) značí kořeny polynomu p, δar (pi ) značí změnu koeficientu ar polynomu pi a Cjr jsou hodnoty absolutních koeficientů podmíněnosti. Srovnejme aproximované hodnoty ξ j (pi ), i = 1, 2, se skutečnými hodnotami kořenů ξj (p1 ) a ξj (p2 ) (viz tabulka 21). V případě polynomu p1 je jeho změna vůči polynomu p natolik malá, že všechny kořeny zůstávají reálné a relativní chyba aproximovaných hodnot kořenů od skutečných hodnot je menší než 10−2 . V případě polynomu p2 je jeho změna vůči polynomu p již natolik závažná, že se dva kořeny změnily v komplexní.
Závěr
50
Závěr V práci je předložena řada numerických metod řešení nelineárních rovnic. Příslušná metoda je vždy podrobně popsána. Je vysloveno tvrzení o konvergenci dané metody, a většinou je toto tvrzení následně dokázáno. Rovněž je uváděno tvrzení (i s důkazem) o řádu konvergence dané metody. Pro lepší názornost byl princip metody vždy prezentován na testovacím příkladu. Každou z metod byl určen tentýž kořen jedné konkrétní nelineární rovnice. Pro demonstraci speciálních vlastností (pomalá konvergence, či nekonvergence) některých metod byly navíc použity i další příklady. Zajímavý je poslední odstavec první kapitoly, který rekapituluje uvedené metody. Metody jsou zde porovnávány z hlediska rychlosti konvergence, a to jednak pomocí řádu konvergence, a jednak pro testovací příklad podle toho, kolik iterací muselo být tou či onou metodou provedeno, aby byl kořen získán s předepsanou přesností. Problém určení všech kořenů polynomu je zde řešen ve dvou etapách. Nejprve se vypočítá největší kořen polynomu. To lze provést např. Newtonovou metodou, nebo lze, jak práce dokládá, s výhodou využít dvojkrokovou metodu. Pro určení ostatních kořenů polynomu jsou předloženy dvě metody. Metoda dělení kořenovým činitelem obecně vede ke zkresleným výsledkům, což je v práci prokázáno a demonstrováno na názorném příkladu. Druhá, Maehlyho metoda, umožňuje určit další kořeny polynomu s libovolnou přesností, a to dokonce i v případě, že dosud určené kořeny nejsou přesné. Poslední kapitola se věnuje studiu citlivosti kořenů polynomu. Z úvah vedených v práci jasně vyplývá, že existují jednak polynomy, jejichž hodnoty kořenů jsou na malou změnu libovolného koeficientu polynomu málo citlivé (tzv. dobře podmíněné polynomy), ale také polynomy, u nichž malá relativní změna koeficientu způsobí mnohem větší relativní změnu v kořenech (tzv. špatně podmíněné polynomy). Je vysloveno a prokázáno tvrzení, že každý polynom, mající kořen násobnosti vyšší než jedna, je špatně podmíněným polynomem. Veškeré výpočty a grafické výstupy byly realizovány v matematickém softwarovém prostředí MATLAB. Takto vzniklo mnoho programů (m–soubory) a obrázků (fig-soubory). Všechny tyto soubory jsou umístěny na CD-ROM, jež je nedílnou součástí rigorózní práce.
Literatura
51
Literatura [1] BULIRSCH, R.; STOER, J. Introduction to Numerical Analysis. New York : Springer-Verlag, 1983. 604 s. ISBN 0-387-90420-4. [2] DATTA, B. N. Numerical Linear Algebra and Applications. International Thomson Publishing Inc., 1995. [3] DĚMIDOVIČ, B. P.; MARON, I. A. Základy numerické matematiky. 1. vyd. Praha : SNTL, 1966. 724 s. [4] DONT, M. Numerické metody – cvičení. Praha : ČVUT, Fakulta elektrotechnická, 1990. 160 s. ISBN 80-01-00400-7. [5] HANSELMAN, D.; LITTLEFIELD, B. Mastering Matlab 6. New Jersey : Prentice Hall, 2001. 832 s. ISBN 0-13-019468-9. [6] HOROVÁ, I. Numerické metody. 1. vyd. Brno : MU v Brně, 1999. 230 s. ISBN 80-210-2202-7. [7] MATHEWS, J. H. Numerical Methods for Mathematics Science and Engineering. New Jersey : Prentice-Hall International, Inc., 1992. [8] ORTEGA, J. M.; RHEINBOLDT, W. C. Iterative Solution of Nonlinear Equations in Several Variables. New York : Academic Press, 1970. [9] PAN, V. Y. Solving a polynomial equation : Some history and recent progress. SIAM Rev., 1997, vol. 39, No. 2, s. 187–220. [10] RALSTON, A. Základy numerické matematiky. 1. vyd. Praha : Academia, 1973. 636 s. [11] ŠÍPEK, J.; ZÍTKO, J. Algoritmy na výpočet kořenů polynomu. Pokroky matematiky, fyziky a astronomie, 2001, roč. 46, č. 1, s. 33–42. [12] ŠMEREK, M. Metody řešení nelineárních rovnic. In XVIII. mezinárodní kolokvium o řízení osvojovacího procesu (sborník příspěvků). Vyškov : VVŠ PV, 2000. s. 337–340. ISBN 80-7231-059-3. [13] ŠMEREK, M. Výpočet kořenů polynomu. In XIX. mezinárodní kolokvium o řízení osvojovacího procesu (sborník příspěvků). Vyškov : VVŠ PV, 2001. s. 399–402. ISBN 80-7231-071-2. [14] ŠMEREK, M. Špatně podmíněné polynomy. In XX. mezinárodní kolokvium o řízení osvojovacího procesu (sborník příspěvků). Vyškov : VVŠ PV, 2002. s. 421–423. ISBN 80-7231-090-9. [15] ŠMEREK, M. Citlivost polynomů a podmíněnost matic. XXI. mezinárodní kolokvium o řízení osvojovacího procesu : sborník abstraktů a elektronických verzí příspěvků na CD-ROMu. [CD-ROM]. Vyškov : VVŠ PV, 2003. ISBN 80-7231-105-0. Adresář: /clanky/1sxmerem.pdf.
Literatura
52
[16] ŠMEREK, M. Citlivost polynomů. XXII. mezinárodní kolokvium o řízení osvojovacího procesu : sborník abstraktů a elektronických verzí příspěvků na CD-ROMu. [CD-ROM]. Vyškov : VVŠ PV, 2004. ISBN 80-7231-116-6. Adresář: /6clanky/1smerekm.pdf. [17] ŠMEREK, M. Citlivost kořenů polynomů. Učitel matematiky. Praha : JČMF, 2005 (přijato k tisku). 11 s. ISSN 1210-9037. [18] VITÁSEK, E. Numerické metody. Praha : SNTL, 1987. [19] WILKINSON, J. H. The evaluation of the zeros of ill-conditioned polynomials. Part I. Numerische Mathematik 1, 1959, vol. 1, s. 150–166. [20] ZAPLATÍLEK, K.; DOŇAR, B. Matlab pro začátečníky. 1. vyd. Praha : BEN, 2003. 144 s. ISBN 80-7300-095-4.
Přílohy
53
Životopis Jmenuji se Michal Šmerek, narodil jsem se 24. 3. 1971 ve Vyškově. Jsem ženatý, mám dceru Michaelu, narozenou 7. 7. 1999. Bydlíme v Zelené Hoře, okres Vyškov. Středoškolské vzdělání jsem získal na gymnáziu Brno, tř. kpt. Jaroše, kde jsem navštěvoval třídu s rozšířenou výukou matematiky. Maturoval jsem dne 25. 5. 1989. V letech 1989–1994 jsem studoval na Přírodovědecké fakultě Masarykovy univerzity v Brně obor Učitelství všeobecně vzdělávacích předmětů matematika – fyzika. Studium jsem ukončil úspěšným složením státní závěrečné zkoušky dne 8. 6. 1994. Dne 1. 10. 1997 jsem zahájil, opět na Přírodovědecké fakultě Masarykovy univerzity v Brně, studium v doktorském studijním programu Obecné otázky matematiky. Od 1. 5. 1998 mi bylo, z rodinných důvodů, studium přerušeno. Po přerušení jsem navázal ve školním roce 1999/2000, opět prvním ročníkem. Státní doktorskou zkoušku jsem složil 29. 5. 2003. Po vykonání roční zakladní vojenské služby v Kroměříži jsem byl, na základě výběrového řízení, dnem 1. září 1995 přijat na místo vysokoškolského učitele Katedry matematiky a informatiky Fakulty ekonomiky obrany státu Vysoké vojenské školy pozemního vojska ve Vyškově. V souvislosti s transformací vysokého vojenského školství byla naše fakulta od 1. září 2005 začleněna do struktur nově vzniklé Univerzity obrany. V současné době tedy působím jako odborný asistent na pozici vedoucího předmětové skupiny Katedry ekonometrie Fakulty ekonomiky a managementu Univerzity obrany v Brně.
V Brně dne 12. 3. 2005
Mgr. Michal Šmerek
Přílohy
54
Publikační činnost Proslovené přednášky: [1] Metody řešení nelineárních rovnic. XVIII. mezinárodní kolokvium o řízení osvojovacího procesu. Vyškov, 18. 5. 2000. [2] Metody řešení nelineárních rovnic. 21. letní škola Historie matematiky. Velké Meziříčí, 28. 8. 2000. [3] Výpočet kořenů polynomu. XIX. mezinárodní kolokvium o řízení osvojovacího procesu. Vyškov, 17. 5. 2001. [4] Výpočet kořenů polynomu. Didaktický seminář. Brno, 6. 5. 2002. [5] Špatně podmíněné polynomy. XX. mezinárodní kolokvium o řízení osvojovacího procesu. Vyškov, 16. 5. 2002. [6] Citlivost polynomů a podmíněnost matic. XXI. mezinárodní kolokvium o řízení osvojovacího procesu. Vyškov, 22. 5. 2003. Publikace: [1] ŠMEREK, M. Metody řešení nelineárních rovnic. In XVIII. mezinárodní kolokvium o řízení osvojovacího procesu (sborník příspěvků). Vyškov : VVŠ PV, 2000. s. 337–340. ISBN 80-7231-059-3. [2] ŠMEREK, M. Výpočet kořenů polynomu. In XIX. mezinárodní kolokvium o řízení osvojovacího procesu (sborník příspěvků). Vyškov : VVŠ PV, 2001. s. 399–402. ISBN 80-7231-071-2. [3] ŠMEREK, M. Špatně podmíněné polynomy. In XX. mezinárodní kolokvium o řízení osvojovacího procesu (sborník příspěvků). Vyškov : VVŠ PV, 2002. s. 421–423. ISBN 80-7231-090-9. [4] ŠMEREK, M. Citlivost polynomů a podmíněnost matic. XXI. mezinárodní kolokvium o řízení osvojovacího procesu : sborník abstraktů a elektronických verzí příspěvků na CD-ROMu. [CD-ROM]. Vyškov : VVŠ PV, 2003. ISBN 80-7231-105-0. Adresář: /clanky/1sxmerem.pdf. [5] ŠMEREK, M. Citlivost polynomů. XXII. mezinárodní kolokvium o řízení osvojovacího procesu : sborník abstraktů a elektronických verzí příspěvků na CD-ROMu. [CD-ROM]. Vyškov : VVŠ PV, 2004. ISBN 80-7231-116-6. Adresář: /6clanky/1smerekm.pdf. [6] ŠMEREK, M. Citlivost kořenů polynomů. Učitel matematiky. Praha : JČMF, 2005 (přijato k tisku). 11 s. ISSN 1210-9037.
Přílohy
55
CD-ROM Přiložený CD-ROM obsahuje: • elektronickou verzi rigorózní práce ve formátu PDF; • m-soubory, vytvořené v matematickém softwarovém prostředí MATLAB. Jsou to programy vytvořené podle v práci uvedených algoritmů. Pomocí nich byly spočteny zde použité příklady. Každý z těchto souborů obsahuje nápovědu, jež uživatele seznámí s daným programem, vysvětlí vstupní a výstupní parametry sledovaného m-souboru. Snadno lze tyto programy použít i pro řešení jiných úloh, některé přímo (změnou parametrů), jiné si mohou vyžádat drobnou úpravu; • fig-soubory. Jsou to obrázky vytvořené opět v matematickém softwarovém prostředí MATLAB. Často jsou vygenerované právě některými z uvedených m-souborů.