Koˇ reny neline´ arn´ıch funkc´ı Matematick´e algoritmy (K611MAG)
Jan Pˇrikryl 9. pˇredn´aˇska 11MAG pondˇel´ı 25. listopadu 2013 verze:2013-11-25 16:47
Obsah 1 Neline´ arn´ı rovnice
1
1.1
Formulace u ´lohy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.2
Existence a jednoznaˇcnost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.3
Podm´ınˇenost ˇreˇsen´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
1.4
Iteraˇcn´ı metody a jejich konvergence . . . . . . . . . . . . . . . . . . . . . . . . .
5
2 Iteraˇ cn´ı metody
7
2.1
Metoda p˚ ulen´ı intervalu neboli bisekce . . . . . . . . . . . . . . . . . . . . . . . .
7
2.2
Metoda postupn´ ych aproximac´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
2.3
Newtonova metoda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.4
Metoda seˇcen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3 Dodatky
16
3.1
Bezpeˇcn´e metody . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2
Numerick´ y v´ ypoˇcet koˇren˚ u polynomu . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.3
Numerick´e ˇreˇsen´ı soustav neline´arn´ıch rovnic . . . . . . . . . . . . . . . . . . . . 17
1 1.1
ˇ sen´ı neline´ Reˇ arn´ıch rovnic Formulace u ´ lohy
Pro detailnˇejˇs´ı obezn´ amen´ı s pojmy, uv´adˇen´ ymi n´ıˇze, doporuˇcuji i zde konzultovat knihu Michaela T. Heathe [3], pˇr´ıpadnˇe nˇejakou z ˇcesk´ ych uˇcebnic ˇci mnoha skript o numerick´e matematice, kter´a v posledn´ıch letech vyˇsla – napˇr´ıklad [1], [2] (ˇc´asti tohoto skripta jsou dostupn´e i on-line).
1
Mnoh´e ze zde pouˇzit´ ych obr´ azk˚ u jsme pˇrevzali pr´avˇe z [3]. Budeme se zde zab´ yvat pˇredevˇs´ım numerick´ ymi metodami pro (pˇribliˇzn´e) ˇreˇsen´ı neline´ arn´ı rovnice f (x) = 0, (1) ˇ sit line´ kde f : R → R je re´ aln´ a neline´ arn´ı funkce jedn´e re´aln´e promˇenn´e. Reˇ arn´ı rovnice tvaru ax + b = c, tj. ax + b − c = 0, jsme se nauˇcili jiˇz na stˇredn´ı ˇskole. Sezn´amili jsme se tam i s ˇreˇsen´ım nˇekter´ ych neline´ arn´ıch rovnic, napˇr´ıklad kvadratick´e rovnice ax2 + bx + c = 0 nebo rovnice sin x = 0. K ˇreˇsen´ı vˇetˇsiny neline´arn´ıch rovnic vˇsak potˇrebujeme pouˇz´ıt nˇekterou vhodnou numerickou metodu. ˇ sen´ım nebo koˇ Reˇ renem uveden´e rovnice nebo tak´e nulov´ ym bodem funkce f naz´ yv´ame takov´e re´aln´e ˇc´ıslo x∗ , pro kter´e plat´ı f (x∗ ) = 0. Neline´arn´ı rovnice mohou m´ıt pr´avˇe jedno ˇreˇsen´ı (rovnice x − sin x = 0), v´ıce ˇreˇsen´ı (rovnice x2 − 1 = 0), nebo nemus´ı m´ıt ˇz´adn´e ˇreˇsen´ı (rovnice sin x = 2). Budeme se tak´e zab´ yvat speci´ aln´ım rovnicemi tvaru x = g(x),
(2)
ˇ sen´ı takov´e rovnice se naz´ kde opˇet g : R → R a tedy f (x) = x − g(x). Reˇ yv´a tak´e pevn´ ym bodem funkce g. Speci´aln´ı situace nast´ av´ a tak´e, je-li uvaˇzovan´a funkce f polynom. Hled´ame pak totiˇz ˇcasto i jeho komplexn´ı koˇreny. Pro polynomy existuj´ı tedy i speci´aln´ı numerick´e metody, jimiˇz se zde vˇsak nem˚ uˇzeme extra zab´ yvat. V praxi se setk´ame i se soustavami neline´arn´ıch rovnic, jejichˇz ˇreˇsen´ım pak nen´ı ˇc´ıslo, ale vektor hodnot. Existuj´ı numerick´e metody i pro ˇreˇsen´ı takov´ ych soustav, z ˇcasov´ ych d˚ uvod˚ u se jimi ale zde tak´e zab´ yvat nebudeme. Pˇr´ıpadn´e z´ajemce odkazujeme na Heathovu knihu [3] nebo na M´ıkovu uˇcebnici [1] ˇci skripta [2].
1.2
Existence a jednoznaˇ cnost
Existence a jednoznaˇcnost ˇreˇsen´ı neline´arn´ıch rovnic je podstatnˇe komplikovanˇejˇs´ı z´aleˇzitost neˇz je tomu u line´ arn´ıch rovnic a jejich soustav. V mnoha pˇr´ıpadech je obt´ıˇzn´e stanovit existenci nebo poˇcet ˇreˇsen´ı neline´ arn´ı rovnice. Zat´ımco u soustav line´arn´ıch rovnic mus´ı b´ yt poˇcet ˇreˇsen´ı roven nule, jedn´e nebo b´ yt nekoneˇcn´ y, neline´arn´ı rovnice mohou m´ıt jak´ ykoli poˇcet ˇreˇsen´ı, a to dokonce i pro jedinou rovnici. Pˇr´ıklad 1 (Existence koˇren˚ u). Uvaˇzujeme-li koˇreny n´asleduj´ıc´ıch rovnic na cel´em R, pak rovnice • ex + 1 = 0 nem´ a ˇz´ adn´e ˇreˇsen´ı • e−x − x = 0 m´ a pr´ avˇe jedno ˇreˇsen´ı • x2 − 4 sin x = 0 m´ a dvˇe ˇreˇsen´ı • x3 + 6x2 + 11x − 6 = 0 m´ a tˇri ˇreˇsen´ı • sin x = 0 m´ a nekoneˇcnˇe mnoho ˇreˇsen´ı Jakkoli je tedy obt´ıˇzn´e z´ıskat jak´ akoli glob´aln´ı tvrzen´ı o poˇctu ˇreˇsen´ı neline´arn´ı rovnice, m´ame pˇresto k dispozici nˇekter´ a uˇziteˇcn´ a lok´aln´ı krit´eria zaruˇcuj´ıc´ı existenci aspoˇ n jednoho koˇrene rovnice na dan´em intervalu. Jedno takov´e praktick´e krit´erium je zaloˇzeno na matematick´e vˇetˇe, kter´a poch´az´ı od B. Bolzana a ˇr´ık´ a: 2
y
x2
x1 a
b x
f(x)
Obr´ azek 1: Pˇr´ıklad neline´ arn´ı funkce s koˇreny x1 a x2 na intervalu [a, b]. Tento interval je pro danou funkci f (x) uz´ avˇerou koˇrene.
Vˇ eta 2 (Bolzanova vˇeta). Necht’ funkce f je spojit´ a na uzavˇren´em omezen´em intervalu [a, b] a necht’ plat´ı f (a)·f (b) < 0 (tj. f (a) a f (b) maj´ı opaˇcn´ a znam´enka). Pak funkce f m´ a na intervalu (a, b) alespoˇ n jeden nulov´y bod. Takov´ y interval [a, b], v jehoˇz koncov´ ych bodech m´a funkce f opaˇcn´a znam´enka, budeme naz´ yvat uz´ avˇ erou ˇreˇsen´ı neline´ arn´ı rovnice f (x) = 0. Jak uvid´ıme pozdˇeji, v ˇradˇe numerick´ ych metod pro ˇreˇsen´ı neline´ arn´ıch rovnic hraje d˚ uleˇzitou roli pr´avˇe postupn´e zuˇzov´an´ı takov´e pˇredem nalezen´e uz´avˇery. Jak ovˇsem poˇc´ ateˇcn´ı uz´avˇeru naj´ıt, je v´ıcem´enˇe z´aleˇzitost´ı pokus˚ u a omyl˚ u. Jedna z moˇznost´ı je odhadnout nˇejak poˇc´ateˇcn´ı interval, na nˇemˇz budeme chov´an´ı funkce f zkoumat (i kdyˇz to jeˇstˇe nebude uz´ avˇera koˇrene), a pak proch´azet t´ımto intervalem po nˇejak´ ych vhodnˇe volen´ ych kroc´ıch, postupnˇe poˇc´ıtat hodnoty f (x) a sledovat, kdy v nich dojde ke zmˇenˇe znam´enka. Poznamenejme jeˇstˇe, ˇze pr´ avˇe zm´ınˇen´e krit´erium ud´av´a pouze postaˇcuj´ıc´ı, nikoli nutnou podm´ınku existence koˇrene. Nemˇelo by n´ as tedy od hled´an´ı koˇrene pˇredem odrazovat to, ˇze jsme nenaˇsli jeho uz´avˇeru. V ˇradˇe pˇr´ıpad˚ u totiˇz ani pro dan´ y koˇren uz´avˇera neexistuje. Jako pˇr´ıklad staˇc´ı vz´ıt trivi´aln´ı rovnici x2 = 0 s jedin´ ym koˇrenem x = 0. Zde pro vˇsechna x m´ame x2 ≥ 0 a ke zmˇenˇe znam´enka tedy nem˚ uˇze doj´ıt. Obrat’me svoji pozornost nyn´ı k rovnici (2). Zde lze k d˚ ukazu existence pevn´eho bodu na dan´em intervalu vyuˇz´ıt opˇet klasick´e matematiky, konkr´etnˇe tzv. vˇety o kontrakci. ˇ Definice 3 (Kontrakce). Rekneme, ˇze funkce g : R → R je na mnoˇzinˇe S ⊆ R kontrakce, pokud existuje konstanta L ∈ R, 0 < L < 1 takov´a, ˇze pro vˇsechna x, y ∈ S plat´ı |g(x) − g(y)| ≤ L|x − y|. Vˇ eta 4 (O existenci pevn´eho bodu). Jestliˇze je funkce g kontrakce na uzavˇren´e mnoˇzinˇe S ⊆ R a g(S) ⊆ S, pak m´ a g v S pevn´y bod, a to pr´ avˇe jeden. Pokud ˇcten´aˇri vad´ı, ˇze jsme v´ yˇse uvedenou definici a vˇetu formulovali pro obecnou mnoˇzinu S, m˚ uˇze si pod S pˇredstavovat nˇejak´ y interval. Z uveden´e vˇety ihned plyne, ˇze pokud v rovnici (1) m˚ uˇzeme ps´ at f (x) = x − g(x), kde g je kontrakce na nˇejak´e uzavˇren´e mnoˇzinˇe S takov´a, ˇze zobrazuje S do sebe sam´e, m´ a rovnice f (x) = 0 na mnoˇzinˇe S pr´avˇe jedno ˇreˇsen´ı, totiˇz pevn´ y 3
bod funkce g. Brzy uvid´ıme, ˇze tato skuteˇcnost n´am d´av´a moˇznost odvodit nˇekter´e numerick´e metody pro ˇreˇsen´ı neline´ arn´ıch rovnic. Poznamenejme jeˇstˇe, ˇze pokud pro vˇsechna x ∈ S existuje derivace funkce g a plat´ı |g 0 (x)| < 1, d´ a se uk´azat, ˇze g na S je kontrakce. Z matematick´eho ˇci logick´eho hlediska jsou naˇse u ´vahy o uz´avˇeˇre spojit´e funkce ˇci pˇredpoklady vˇety o kontrakci pouze postaˇcuj´ıc´ı podm´ınky, nikoli vˇsak podm´ınky nutn´e. Nen´ı tedy nikde ps´ano, ˇze funkce f nem˚ uˇze m´ıt nulov´ y bod v intervalu, kter´ y nen´ı uz´avˇerou, nebo ˇze funkce g, kter´a nen´ı kontrakce, nem˚ uˇze m´ıt pevn´ y bod. V praxi tedy m˚ uˇze b´ yt uˇziteˇcn´e vyuˇz´ıt soudob´ ych moˇznost´ı naˇseho softwarov´eho vybaven´ı a pˇri ˇreˇsen´ı neline´arn´ıch rovnic si nejprve nechat vykreslit graf funkce f pro rovnici (1) nebo grafy y = x a y = g(x) pro rovnici (2). • funkce f (x) = x2 m´a nulov´ y bod x = 0 na
Pˇr´ıklad 5 (Nesplnˇen´e postaˇcuj´ıc´ı podm´ınky). intervalu [−1, 1], kter´ y nen´ı uz´ avˇera
• funkce g(x) = sin x m´ a pevn´ y bod x = 0, ale v okol´ı tohoto bodu to nen´ı kontrakce Doposud jsme se soustˇredili pˇrev´ aˇznˇe na existenci koˇren˚ u neline´arn´ıch rovnic a ne na jejich jednoznaˇcnost, protoˇze se obecnˇe m´ a za to, ˇze neline´arn´ı rovnice mohou m´ıt v´ıce neˇz jedno ˇreˇsen´ı, pˇrinejmenˇs´ım glob´ alnˇe. Jednoznaˇcnost koˇrene n´as pˇresto m˚ uˇze zaj´ımat, alespoˇ n lok´alnˇe, napˇr´ıklad na dan´em intervalu. Pˇripomeˇ nme si, ˇze z line´arn´ı algebry v´ıme, ˇze soustava line´arn´ıch rovnic s regul´ arn´ı matic´ı m´ a vˇzdy pr´ avˇe jedno ˇreˇsen´ı. Pro neline´arn´ı funkce f plat´ı podobn´e tvrzen´ı o regularitˇe, pˇrinejmenˇs´ım lok´ alnˇe. Pokud totiˇz funkce f m´a v dan´em bodˇe x∗ nenulovou derivaci, pak existuje otevˇren´ y interval kolem tohoto bodu, v nˇemˇz je funkce f ostˇre monot´onn´ı, tedy rostouc´ı nebo klesaj´ıc´ı. V takov´e situaci v okol´ı bodu x∗ tedy m˚ uˇze existovat nejv´ yˇse jeden koˇren. Pokud vˇsak v nˇejak´em koˇrenu x∗ m´ a funkce f nulovou derivaci, m´a tento koˇren jist´e zvl´aˇstn´ı vlastnosti, kter´e ovlivˇ nuj´ı jak podm´ınˇenost ˇreˇsen´e u ´lohy, tak tak´e chov´an´ı pouˇzit´e numerick´e metody. Nulov´ y bod x∗ funkce f , pro kter´ y plat´ı z´aroveˇ n f (x∗ ) = 0 a f 0 (x∗ ) = 0 se naz´ yv´ a n´ asobn´y koˇren rovnice f (x) = 0. Geometricky to znamen´a, ˇze graf funkce f m´a v tomto bodˇe vodorovnou teˇcnu spl´ yvaj´ıc´ı s osou x. Koˇreny, kter´e nejsou n´asobn´e, se naz´ yvaj´ı jednoduch´e. Koˇren x1 z Obr´ azku1 je tedy jednoduch´ y, koˇren x2 je n´asobn´ y. Pojem n´asobnosti lze pro hladk´e funkce f d´ale upˇresnit. Pokud plat´ı f (x∗ ) = f 0 (x∗ ) = f 00 (x∗ ) = · · · = f (m−1) (x∗ ) = 0, ale f (m) 6= 0, ˇrekneme, ˇze n´ asobnost koˇrene x∗ je m. Pˇr´ıklad 6 (Pevn´e body). Ovˇeˇrte si n´ asleduj´ıc´ı tvrzen´ı: • funkce g(x) = sin x m´ a jedin´ y pevn´ y bod x = 0 • funkce g(x) = x2 m´ a pevn´e body x = 0 a x = 1 • koˇren x = 0 rovnice sin x = 0 je jednoduch´ y • koˇren x = 0 rovnice x2 = 0 je dvojn´asobn´ y • koˇren x = 1 rovnice x3 − 3x2 + 3x − 1 = 0 je trojn´asobn´ y
1.3
Podm´ınˇ enost ˇ reˇ sen´ı
Abychom mohli kvantitativnˇe mˇeˇrit citlivost ˇreˇsen´ı neline´arn´ıch rovnic na data (funkˇcn´ı hodnoty), mus´ıme pracovat s absolutn´ım ˇc´ıslem podm´ınˇenosti, coˇz je obdoba jiˇz zaveden´eho ˇc´ısla podm´ınˇenosti z minul´e pˇredn´ aˇsky, kde ale m´ısto relativn´ıch zmˇen v ˇcitateli i jmenovateli vystupuj´ı absolutn´ı odchylky. Je to d´ ano t´ım, ˇze hodnota funkce f v koˇrenu rovnice je rovna nule. 4
Obr´ azek 2: Podm´ınˇenost koˇren˚ u neline´ arn´ı rovnice f (x) = 0. Vlevo: dobˇre podm´ınˇen´a u ´loha, vpravo: ˇspatnˇe podm´ınˇen´ au ´loha.
D´a se uk´azat, ˇze pokud m´ a funkce f v okol´ı koˇrene x∗ derivaci, je pak toto ˇc´ıslo podm´ınˇenosti pˇribliˇznˇe 1 Cp,abs ≈ 0 ∗ . |f (x )| Pokud je ve jmenovateli f 0 (x∗ ) = 0 (n´asobn´ y koˇren), klademe Cp,abs = ∞. Z definice ˇc´ısla podm´ınˇenosti pak vypl´ yv´ a, ˇze pokud najdeme bod x ˜ takov´ y, ˇze |f (˜ x)| < , m˚ uˇze odchylka |˜ x − x∗ | tohoto bodu od koˇrene rovnice f (x) = 0 m´ıt velikost /|f 0 (x∗ )|. Pro mal´e hodnoty |f 0 (x∗ )| tedy m˚ uˇze b´ yt tato odchylka od koˇrene velk´a, i kdyˇz funkˇcn´ı hodnota sama je mal´a. ˇ arkovan´e kˇrivky vyznaˇcuj´ı oblast nejistoty kolem kaˇzd´e plnˇe Celou situaci ilustruje Obr´ azek 2. C´ nakreslen´e kˇrivky, takˇze nulov´ y bod dan´e funkce m˚ uˇze b´ yt kdekoli mezi body, v nichˇz ˇc´arkovan´e kˇrivky prot´ınaj´ı vodorovnou osu. Mal´ y interval nejistoty pro nulov´ y bod na lev´em obr´azku je d´an t´ım, ˇze dan´ a kˇrivka strmˇe roste (takˇze pˇrevr´acen´a hodnota derivace je mal´a), kdeˇzto velk´ y interval nejistoty pro nulov´ y bod na prav´em obr´azku plyne z pomal´eho r˚ ustu (a tedy velk´e pˇrevr´acen´e hodnoty derivace). Vˇsimnˇete si tak´e toho, ˇze ˇs´ıˇre p´asu nejistot kolem funkˇcn´ıch hodnot je na obou obr´ azc´ıch stejn´ a. M´ame-li n´asobn´ y koˇren x∗ , je f 0 (x∗ ) = 0, takˇze ˇc´ıslo podm´ınˇenosti n´asobn´eho koˇrene je nekoneˇcn´e. To d´ av´ a smysl, protoˇze nepatrn´a zmˇena v f m˚ uˇze zp˚ usobit, ˇze z n´asobn´eho koˇrene se stane v´ıce neˇz jeden koˇren nebo naopak n´asobn´ y koˇren zmiz´ı. Staˇc´ı si k tomu nakreslit napˇr´ıklad funkci f (x) = x2 a posunout ji o mal´e nahoru nebo dol˚ u. pˇr´ıklad,obr´ azek Podm´ınˇenost neline´ arn´ı rovnice ovlivˇ nuje n´aˇs pohled na pˇribliˇzn´e ˇreˇsen´ı x ˜: m´ame usilovat o to, aby |f (˜ x)| byla mal´ a, nebo sp´ıˇse o to, aby bylo mal´e |˜ x − x∗ |, jakkoli pˇresn´e ˇreˇsen´ı x∗ pˇredem nezn´ ame? Jak uˇz to u numerick´ ych metod b´ yv´a, obˇe uveden´e veliˇciny nejsou nutnˇe mal´e souˇcasnˇe, z´ avis´ı to jeˇstˇe na podm´ınˇenosti. Tato skuteˇcnost ovlivˇ nuje volbu algoritm˚ u numerick´ ych metod, o nichˇz budeme hovoˇrit ve zbytku t´eto pˇredn´aˇsky. V kaˇzd´em pˇr´ıpadˇe je uˇziteˇcn´e z´ıskat pˇredem nˇejakou informaci o podm´ınˇenosti ˇreˇsen´e u ´lohy.
1.4
Iteraˇ cn´ı metody a jejich konvergence
Numerick´e metody pro ˇreˇsen´ı neline´ arn´ıch rovnic jsou vesmˇes metody iteraˇcn´ı. Iterace (z lat. iterare, opakovat) znamen´ a postupn´e opakov´an´ı urˇcit´eho postupu, bˇehem kter´eho se postupnˇe generuje posloupnost hodnot x0 , x1 , . . . , xk , . . . takov´a, ˇze v naˇsem pˇr´ıpadˇe (hled´ame koˇren x∗ neline´arn´ı rovnice) postupnˇe z´ısk´ avan´e hodnoty konverguj´ı k hledan´emu ˇreˇsen´ı, xk → x∗ pro k → ∞. Pˇri skuteˇcn´em v´ ypoˇctu samozˇrejmˇe nem˚ uˇzeme j´ıt s k do nekoneˇcna a iteraˇcn´ı postup zastav´ıme po urˇcit´em dostateˇcnˇe velk´em poˇctu krok˚ u pomoc´ı vhodnˇe zvolen´eho zastavovac´ıho
5
krit´eria. Z´ısk´ ame tak pˇribliˇznou hodnotu hledan´eho koˇrene. Term´ın iterace se v numerick´e matematice pouˇz´ıv´ a nejen k oznaˇcen´ı v´ yˇse uveden´eho postupu jako celku, ale naz´ yv´a se tak tak´e kaˇzd´ y jeho krok a naz´ yvaj´ı se tak tak´e postupnˇe poˇc´ıtan´e hodnoty xk , tedy aproximace hledan´eho koˇrene. Abychom mohli porovn´ avat efektivitu iteraˇcn´ıch metod, potˇrebujeme nˇejak charakterizovat jejich rychlost konvergence, tj. rychlost konvergence posloupnosti iterac´ı xk k hledan´emu koˇrenu rovnice. Chyba(nepˇresnost) k-t´e iterace, kterou budeme oznaˇcovat ek , se obvykle definuje jako ek = xk − x∗ , kde xk je aproximace (pˇribl´ıˇzen´ı) hledan´eho ˇreˇsen´ı z´ıskan´a v iteraci k a x∗ je skuteˇcn´e (pˇresn´e) ˇreˇsen´ı. Nˇekter´e z pouˇz´ıvan´ ych metod neprodukuj´ı pˇr´ımo konkr´etn´ı pˇribliˇzn´e ˇreˇsen´ı xk , ale pouze interval, kter´ y s urˇcitost´ı obsahuje pˇresn´e ˇreˇsen´ı, pˇriˇcemˇz d´elka tohoto intervalu se bˇehem iteraˇcn´ıho procesu postupnˇe zmenˇsuje. U takov´e metody pak ek definujeme jako d´elku tohoto intervalu po k-t´e iteraci. V obou pˇr´ıpadech pak ˇrekneme, ˇze dan´a iteraˇcn´ı metoda konverguje s rychlost´ı r (tak´e: metoda je ˇr´adu r), jestliˇze pro nˇejakou koneˇcnou kladnou konstantu C > 0 plat´ı |ek+1 | lim = C. k→∞ |ek |r Speci´alnˇe se rozliˇsuj´ı n´ asleduj´ıc´ı pˇr´ıpady: • pokud r = 1 a C < 1, je konvergence line´ arn´ı, • pokud r > 1, je konvergence superline´ arn´ı, • pokud r = 2, je konvergence kvadratick´ a, • pokud r = 3, je konvergence kubick´ a, atd. Jeden z d˚ uvod˚ u, proˇc rozliˇsujeme mezi line´arn´ı a superline´arn´ı konvergenc´ı, je ten, ˇze, asymptoticky pro velk´ a k, line´ arnˇe konvergentn´ı posloupnost z´ısk´av´a po kaˇzd´e iteraci jist´ y st´ale stejn´ y poˇcet pˇresn´ ych ˇc´ıslic, kdeˇzto superline´arnˇe konverguj´ıc´ı posloupnost v jednotliv´ ych iterac´ıch z´ısk´av´a poˇcet pˇresn´ ych ˇc´ıslic, kter´ y st´ ale roste. Pˇresnˇeji m˚ uˇzeme ˇr´ıci, ˇze line´arnˇe konvergentn´ı posloupnost z´ısk´ av´ a v kaˇzd´e iteraci − log(C) pˇresn´ ych ˇc´ıslic, kdeˇzto superline´arnˇe konverguj´ıc´ı posloupnost m´ a po kaˇzd´e iteraci r-kr´ at v´ıce pˇresn´ ych ˇc´ıslic neˇz mˇela pˇred touto iterac´ı. Speci´alnˇe pak plat´ı, ˇze u kvadraticky konvergentn´ı metody se poˇcet pˇresn´ ych ˇc´ıslic po kaˇzd´e iteraci zdvojn´asob´ı (pro dostateˇcnˇe velk´ a k). Pˇr´ıklad 7 (Rychlosti konvergence). Jestliˇze ˇcleny n´asleduj´ıc´ıch posloupnost´ı pˇredstavuj´ı velikosti chyb postupnˇe generovan´ ych iteraˇcn´ıch aproximac´ı, jsou rychlosti konvergence takov´e, jak je u jednotliv´ ych posloupnost´ı uvedeno. • 10−2 , 10−3 , 10−4 , 10−5 , . . .
line´ arn´ı, C = 10−1
• 10−2 , 10−4 , 10−6 , 10−8 , . . .
line´ arn´ı, C = 10−2
• 10−2 , 10−3 , 10−5 , 10−8 , . . .
superline´arn´ı, ale ne kvadratick´a
• 10−2 , 10−4 , 10−8 , 10−16 , . . .
kvadratick´a
V teorii numerick´ ych metod se dokazuj´ı vˇety o konvergenci, kter´e n´am umoˇzn ˇuj´ı ˇr´ıci, kdy pro danou rovnici ta ˇci ona metoda konverguje a jak rychle. Ned´avaj´ı n´am ale explicitnˇe pokyny pro to, kdy m´ ame iteraˇcn´ı proces zastavit a prohl´asit v´ ysledn´e pˇribliˇzn´e ˇreˇsen´ı za dostateˇcnˇe ” pˇresn´e“. Navrhnout vhodn´e zastavovac´ı krit´erium je pomˇernˇe sloˇzit´a z´aleˇzitost, a to z ˇrady 6
d˚ uvod˚ u. D´ıky teorii m˚ uˇzeme v z´ asadˇe vˇedˇet, ˇze se chyba |ek | postupnˇe zmenˇsuje, ale protoˇze nezn´ame pˇresn´e ˇreˇsen´ı, nen´ı tu moˇznost pˇr´ımo zjistit, jak velik´e |ek | je. Rozumnou n´ahraˇzkou tu m˚ uˇze slouˇzit relativn´ı zmˇena v postupn´ ych iterac´ıch, tedy |xk+1 − xk | . |xk | Pokud se tato veliˇcina stane dostateˇcnˇe malou, znamen´a to, ˇze se pˇribliˇzn´e hodnoty ˇreˇsen´ı uˇz pˇrestaly v´ yznamnˇe mˇenit a nem´ a tedy cenu pokraˇcovat. Na druh´e stranˇe bychom chtˇeli m´ıt jistotu, ˇze jsme skuteˇcnˇe z´ıskali dobr´e pˇribliˇzn´e ˇreˇsen´ı a ˇze tedy aspoˇ n je hodnota f (xk ) pˇrimˇeˇrenˇe mal´a. Jak uˇz jsme si ale mohli povˇsimnout, obˇe dvˇe tyto uveden´e veliˇciny nemus´ı b´ yt mal´e souˇcasnˇe, roli tu hraje podm´ınˇenost u ´lohy. D´ale se zde projevuje tak´e pˇr´ıpadn´a zmˇena mˇeˇr´ıtka u promˇenn´e x a funkce f . Ze vˇsech tˇechto d˚ uvod˚ u je vytvoˇren´ı zcela spolehliv´eho zastavovac´ıho krit´eria velmi obt´ıˇzn´e a mus´ıme se tak´e spol´ehat na dalˇs´ı informace, kter´e o ˇreˇsen´e u ´loze v´ıme. U iteraˇcn´ıch metod, kter´e budeme vz´apˇet´ı v t´eto pˇredn´aˇsce popisovat, proto zpravidla vynech´av´ame jak´ ykoli test na konvergenci a m´ısto toho pouze naznaˇcujeme jist´ y neurˇcen´ y poˇcet iterac´ı s t´ım, ˇze iteraˇcn´ı proces je tˇreba ukonˇcit pot´e, co se vyhov´ı urˇcit´emu vhodn´emu krit´eriu, jehoˇz volba je (bohuˇzel) na uˇzivateli.
2
Numerick´ e metody ˇ reˇ sen´ı neline´ arn´ıch rovnic
Budeme se zab´ yvat numerick´ ym (pˇribliˇzn´ ym) ˇreˇsen´ım neline´arn´ı rovnice (1): pro danou spojitou ∗ funkci f : R → R hled´ ame bod x ∈ R takov´ y, ˇze f (x∗ ) = 0.
2.1
Metoda p˚ ulen´ı intervalu neboli bisekce
V poˇc´ıtaˇcov´e aritmetice s koneˇcnou pˇresnost´ı nemus´ı existovat strojov´e ˇc´ıslo x∗ takov´e, ˇze f (x∗ ) je pˇresnˇe nula. Alternativn´ı moˇznost je hledat nˇejak´ y velmi mal´ y interval [a, b], ve kter´em f mˇen´ı znam´enko. Jak jsme jiˇz uvedli v odstavci 1.2, takov´a uz´ avˇera zaruˇcuje, ˇze pˇr´ısluˇsn´a spojit´a funkce mus´ı nˇekde uvnitˇr tohoto intervalu m´ıt nulov´ y bod. Metoda p˚ ulen´ı intervalu neboli metoda bisekce zaˇc´ın´ a od nˇejak´e poˇc´ ateˇcn´ı uz´avˇery a postupnˇe sniˇzuje jej´ı velikost do t´e doby, aˇz je ˇreˇsen´ı uzavˇreno s poˇzadovanou pˇresnost´ı (resp. tak, jak to aritmetika poˇc´ıtaˇce dovol´ı). V kaˇzd´e iteraci se nejprve stanov´ı stˇred aktu´ aln´ıho intervalu a pro dalˇs´ı iteraci se ponech´a pouze jedna z polovin intervalu podle toho, jak´e znam´enko m´a funkˇcn´ı hodnota ve stˇredu. Tato polovina pak tvoˇr´ı opˇet (jiˇz kratˇs´ı) uz´ avˇer, s n´ımˇz vstupujeme do dalˇs´ı iterace. Metodu bisekce form´alnˇe m˚ uˇzeme zapsat jako Algoritmus 1, v nˇemˇz jako vstupn´ı data figuruje funkce f , uz´avˇera [a, b] a chybov´a tolerance ∆tol pro d´elku v´ ysledn´eho intervalu obsahuj´ıc´ıho koˇren. Uvedeme jeˇstˇe p´ ar pozn´ amek k implementaci metody bisekce ve v´ yˇse uveden´em algoritmu: Pˇredevˇs´ım, zd´ a se, ˇze nejpˇrirozenˇejˇs´ı vzorec pro v´ ypoˇcet stˇredu intervalu [a, b] by byl m = (a+b)/2. Jenˇze v poˇc´ıtaˇcov´e aritmetice nen´ı v extr´emn´ıch pˇr´ıpadech zaruˇceno, ˇze takto poˇc´ıtan´ y bod m v˚ ubec padne do intervalu [a, b]. Komu se to zd´a divn´e, m˚ uˇze si v aritmetice se dvˇema des´ıtkov´ ymi ˇc´ıslicemi zkusit podle tohoto vzorce spoˇc´ıtat stˇred intervalu [0.67, 0.69] (vyjde m = 0.7). Kromˇe toho m˚ uˇze u tohoto vzorce meziv´ ysledek a + b pˇrekroˇcit rozsah poˇc´ıtaˇce i v situac´ıch, kdy stˇred intervalu v rozsahu poˇc´ıtaˇce leˇz´ı. Jakkoli jde o extr´emn´ı pˇr´ıpady, je na tomto jednoduch´em pˇr´ıkladu vidˇet, ˇze poˇc´ıtaˇcov´a implementace algoritm˚ u nen´ı jen pouh´e pˇrepisov´an´ı vzoreˇck˚ u v nˇejak´em vhodn´em programovac´ım jazyce. Vzorec pouˇzit´ y v Algoritmu 1 se uveden´ ym probl´em˚ um vyh´ yb´ a. 7
Algoritmus 1 Metoda p˚ ulen´ı intervalu Require: Funkce f , uz´ avˇera [a, b], chybov´a tolerance ∆tol Ensure: x∗ : f (x∗ ) ≈ 0 while (b − a) > ∆tol do b−a m←a+ 2 if sgn f (a) = sgn f (m) then a←m else b←m end if end while D´ale, pokud jde o testov´ an´ı toho, zda dvˇe hodnoty f (x1 ) a f (x2 ) maj´ı stejn´e znam´enko, je na poˇc´ıtaˇci lepˇs´ı pouˇz´ıvat funkci signum neˇz matematicky ekvivalentn´ı testov´an´ı, zda souˇcin f (x1 ) · f (x2 ) je kladn´ y nebo z´ aporn´ y. Takov´ y souˇcin m˚ uˇze totiˇz tak´e pˇrekroˇcit rozsah poˇc´ıtaˇce smˇerem k nekoneˇcnu a v okol´ı koˇrene smˇerem k nule. Poznamenejme pro poˇr´adek, ˇze je sgn(x) = 1 pro x ≥ 0 a sign(x) = −1 pro x < 0. Pˇr´ıklad 8 (Metoda bisekce). Metodu p˚ ulen´ı intervalu uk´aˇzeme na pˇr´ıkladu hled´an´ı koˇrene rovnice f (x) = x2 − 4 sin x = 0. Jako poˇc´ateˇcn´ı uz´ avˇeru vezmeme interval [a, b], kde a = 1 a b = 3. Z´aleˇz´ı tu pouze na tom, aby se funkˇcn´ı hodnoty v tˇechto dvou bodech liˇsily ve znam´enku. Vypoˇc´ıt´ame hodnotu funkce ve stˇredn´ım bodˇe intervalu, tedy v m = 2 a zjist´ıme, ˇze f (m) m´a opaˇcn´e znam´enko neˇz f (a), takˇze si podrˇz´ıme levou polovinu poˇc´ ateˇcn´ıho intervalu a poloˇz´ıme pro dalˇs´ı krok b = m. Pak tento postup opakujeme tak dlouho, aˇz se interval uz´avˇery z´ uˇz´ı na poˇzadovanou velikost. N´asleduj´ıc´ı tabulka ukazuje moˇznou posloupnost iterac´ı. a
f (a)
b
f (b)
1.000000 1.000000 1.500000 1.750000 1.875000 1.875000 1.906250 1.921875 1.929688 1.933594 1.933594 1.933594 1.933594
-2.365884 -2.365884 -1.739980 -0.873444 -0.300718 -0.300718 -0.143255 -0.062406 -0.021454 -0.000846 -0.000846 -0.000846 -0.000846
3.000000 2.000000 2.000000 2.000000 2.000000 1.937500 1.937500 1.937500 1.937500 1.937500 1.935547 1.934570 1.934082
8.435520 0.362810 0.362810 0.362810 0.362810 0.019849 0.019849 0.019849 0.019849 0.019849 0.009491 0.004320 0.001736
Interval, u kter´eho jsme iterace ukonˇcili, m´a d´elku menˇs´ı neˇz 0.0005 a m˚ uˇzeme tedy ˇr´ıci, ˇze nalezen´ y koˇren je s touto pˇresnost´ı roven x∗ ≈ 1.934. Na z´avˇer jeˇstˇe nˇekolik pozn´ amek k metodˇe p˚ ulen´ı intervalu. 8
• V metodˇe p˚ ulen´ı se nikde nevyuˇz´ıvaj´ı velikosti funkˇcn´ıch hodnot, pouze jejich znam´enka. • Pokud v´ ypoˇcet zaˇcneme s uz´ avˇerou spojit´e funkce, pak metoda konverguje vˇzdy, ale dosti pomalu. • V kaˇzd´e iteraci se d´elka uz´ avˇery sniˇzuje na polovinu, takˇze rychlost konvergence je line´ arn´ı, s r = 1 a C = 0.5. • V kaˇzd´e iteraci bisekce z´ısk´ av´ ame jednu dalˇs´ı pˇresnou dvojkovou ˇc´ıslici v pˇribliˇzn´em ˇreˇsen´ı. • Pro dan´ y poˇc´ ateˇcn´ı interval [a, b] je d´elka intervalu po k iterac´ıch rovna (b − a)/2k , takˇze k dosaˇzen´ı chybov´e tolerance tol je zapotˇreb´ı zhruba
log2
b−a tol
iterac´ı, nez´ avisle na vlastnostech pouˇzit´e funkce f .
2.2
Metoda postupn´ ych aproximac´ı
Metoda postupn´ych aproximac´ı nebo tak´e metoda prost´e iterace slouˇz´ı k hled´an´ı pevn´ ych bod˚ u funkce g z rovnice (2). Pˇripomeˇ nme tedy, ˇze pro funkci g : R → R se pevn´ym bodem naz´ yv´ a takov´e ˇc´ıslo x∗ (pokud existuje), pro kter´e plat´ı x∗ = g(x∗ ). D˚ uvodem tohoto n´ azvu je skuteˇcnost, ˇze x∗ se po aplikaci funkce g nezmˇen´ı. Zat´ımco u neline´arn´ı rovnice f (x) = 0 hled´ ame bod, v nˇemˇz graf funkce f prot´ın´a osu x (tedy pˇr´ımku y = 0), pˇri hled´an´ı pevn´eho bodu funkce g chceme naj´ıt bod, v nˇemˇz graf funkce g protne diagon´aln´ı pˇr´ımku ´ y = x. Ulohy na hled´ an´ı pevn´eho bodu dost ˇcasto poch´azej´ı pˇr´ımo z praxe, ale pro n´as zde maj´ı v´ yznam tak´e z toho d˚ uvodu, ˇze ˇreˇsen´ı neline´arn´ı rovnice (1) lze zpravidla pˇrev´est na hled´an´ı pevn´eho bodu odpov´ıdaj´ıc´ı neline´ arn´ı funkce g, tedy na ˇreˇsen´ı rovnice (2). Metoda postupn´ ych aproximac´ı (prost´e iterace) pro ˇreˇsen´ı t´eto rovnice je zaloˇzena na opakovan´em (iteraˇcn´ım) pouˇzit´ı vzorce xk+1 = g(xk ) s vhodnˇe zvolen´ ym poˇc´ ateˇcn´ım pˇribl´ıˇzen´ım (poˇc´ateˇcn´ı aproximac´ı) x0 . Chceme-li ˇreˇsit rovnici f (x) = 0 metodou postupn´ ych aproximac´ı, pak ji nejprve mus´ıme pˇrev´est na u ´lohu o pevn´em bodu pro nˇejakou vhodnˇe vybranou funkci g. Takov´ ych moˇznost´ı b´ yv´ a pro danou f v´ıce, ale ne vˇsechny jsou stejnˇe vhodn´e pro z´ısk´an´ı iteraˇcn´ıho sch´ematu k ˇreˇsen´ı v´ ychoz´ı rovnice. V´ ysledn´ a iteraˇcn´ı metoda se pro r˚ uzn´e volby g m˚ uˇze liˇsit nejen co do rychlosti konvergence, ale tak´e v tom, zda v˚ ubec konverguje ˇci nikoli. ´ Pˇr´ıklad 9 (Ulohy na pevn´ y bod). Neline´arn´ı rovnice f (x) = x2 − x − 2 = 0 m´a koˇreny x∗ = 2 a x∗ = −1. Mezi ekvivalentn´ı u ´lohy na hled´an´ı pevn´eho bodu patˇr´ı u ´lohy (2) s funkcemi (ovˇeˇrte si to) 1. g(x) = x2 − 2, √ 2. g(x) = x + 2 (ekvivalence pouze pro nez´aporn´e pevn´e body, srv. (2)), 9
Obr´ azek 3: Pevn´ y bod (2, 2) neline´arn´ıch funkc´ı.
Obr´ azek 4: Metoda postupn´ ych aproximac´ı pro prvn´ı a druhou funkci g.
3. g(x) = 1 + (2/x), 4. g(x) = (x2 + 2)/(2x − 1). Na obr. 3 je vykreslen pr˚ ubˇeh kaˇzd´e z tˇechto funkc´ı spolu s pˇr´ımkou y = x. Vˇsimnˇeme si, ˇze funkce g jsou konstruov´ any tak, ˇze jejich grafy vesmˇes prot´ınaj´ı pˇr´ımku y = x v pevn´em bodˇe (2, 2). Pr˚ ubˇeh pˇr´ısluˇsn´ ych iteraˇcn´ıch sch´emat metody postupn´ ych aproximac´ı je graficky zn´azornˇen na ˇ Obr´azc´ıch 4 a 5. Sipka ve svisl´em smˇeru odpov´ıd´a v´ ypoˇctu hodnoty dan´e funkce v nˇejak´em bodˇe a vodorovn´a ˇsipka smˇeˇruj´ıc´ı k pˇr´ımce y = x vyznaˇcuje, ˇze se v´ ysledek pˇredchoz´ıho v´ ypoˇctu hodnoty funkce g pouˇzije jeko argument pro pˇr´ıˇst´ı v´ ypoˇcet funkˇcn´ı hodnoty. U prvn´ı z uveden´ ych funkc´ı vid´ıme, ˇze i pˇres to, ˇze poˇc´ ateˇcn´ı bod je velmi bl´ızko ˇreˇsen´ı, postupn´e aproximace diverguj´ı. U ostatn´ıch tˇr´ı funkc´ı je vidˇet, ˇze postupn´e iterace konverguj´ı k pevn´emu bodu, i kdyˇz byly odstartov´any v bodˇe, kter´ y je od ˇreˇsen´ı relativnˇe daleko. Zd´a se pˇritom, ˇze rychlosti konvergence pro tyto tˇri funkce se mohou liˇsit. 10
Obr´ azek 5: Metoda postupn´ ych aproximac´ı pro tˇret´ı a ˇctvrtou funkci g.
Jak lze z graf˚ u funkc´ı na Obr´ azc´ıch 4 a 5 vidˇet, chov´an´ı metody prost´e iterace se m˚ uˇze znaˇcnˇe odliˇsovat, od divergence pˇres pomalou konvergenci k rychl´e konvergenci. Nejjednoduˇsˇs´ı (i kdyˇz ne nejobecnˇejˇs´ı) zp˚ usob, jak charakterizovat chov´an´ı iteraˇcn´ıho sch´ematu xk+1 = g(xk ) pro ˇreˇsen´ı u ´lohy na pevn´ y bod tvaru x = g(x), je pokusit se vz´ıt v u ´vahu derivaci funkce g v hledan´em ˇreˇsen´ı x∗ za pˇredpokladu, ˇze funkce g je hladk´a a tato derivace existuje. D´a se uk´azat, ˇze pokud x∗ = g(x∗ ) a |g 0 (x∗ )| < 1, pak iteraˇcn´ı sch´ema metody postupn´ ych aproximac´ı lok´ alnˇ e ∗ konverguje. To znamen´ a, ˇze existuje nˇejak´ y interval obsahuj´ıc´ı x takov´ y, ˇze metoda prost´e iterace s funkc´ı g konverguje, pokud je odstartov´ana z nˇejak´eho x0 , jeˇz leˇz´ı uvnitˇr tohoto inˇ ık´ tervalu. R´ ame tak´e, ˇze metoda konverguje pro dostateˇcnˇe bl´ızk´e poˇc´ateˇcn´ı pˇribl´ıˇzen´ı. Naproti tomu pokud |g 0 (x∗ )| > 1, pak metoda prost´e iterace diverguje pro jak´ekoli poˇc´ateˇcn´ı pˇribl´ıˇzen´ı kromˇe x∗ . D˚ ukaz tohoto tvrzen´ı je zaloˇzen na vˇetˇe o stˇredn´ı hodnotˇe funkce, ale z ˇcasov´ ych d˚ uvod˚ u jej zde neuv´ad´ıme, jakkoli nen´ı sloˇzit´ y (viz [1], [2] nebo [3]). Plyne z nˇej ale tak´e to, ˇze pokud metoda ˇ ım konverguje, je jej´ı asymptotick´ a rychlost konvergence line´arn´ı s konstantou C = |g 0 (x∗ )|. C´ menˇs´ı je tato konstanta, t´ım je konvergence rychlejˇs´ı, a ide´aln´ı by tedy pro danou rovnici (1) bylo naj´ıt ekvivalentn´ı formulaci (2) s funkc´ı g, pro niˇz by platilo g 0 (x∗ ) = 0. V takov´em pˇr´ıpadˇe se d´ a pomoc´ı Taylorova rozvoje opˇet pomˇernˇe snadno uk´azat, ˇze konvergence je nejm´enˇe kvadratick´ a. V pˇr´ıˇst´ım odstavci si pop´ıˇseme jeden systematick´ y zp˚ usob takov´e volby funkce g pro rovnici f (x) = 0. Pˇr´ıklad 10 (Konvergence metody postupn´ ych aproximac´ı). Pro ˇctyˇri u ´lohy na pevn´ y bod z pˇredch´azej´ıc´ıho pˇr´ıkladu m´ ame n´ asleduj´ıc´ı v´ ysledky: 1. g 0 (x) = 2x, takˇze g 0 (2) = 4 a metoda postupn´ ych aproximac´ı tedy diverguje. √ 2. g 0 (x) = 1/(2 x + 2), takˇze g 0 (2) = 1/4 a metoda postupn´ ych aproximac´ı konverguje lin´arnˇe s konstantou C = 1/4. Kladn´e znam´enko derivace g 0 (2) vede k tomu, ˇze se iterace pˇribliˇzuj´ı k pevn´emu bodu z jedn´e strany. 3. g 0 (x) = −2/x2 , takˇze g 0 (2) = −1/2 a metoda postupn´ ych aproximac´ı konverguje line´arnˇe s konstantou C = 1/2. Z´ aporn´e znam´enko derivace g 0 (2) vede k tomu, ˇze se iterace pˇribliˇzuj´ı k pevn´emu bodu po spir´ ale, stˇr´ıdavˇe vˇzdy z opaˇcn´e strany. 4. g 0 (x) = (2x2 − 2x − 4)/(2x − 1)2 , takˇze g 0 (2) = 0 a metoda postupn´ ych aproximac´ı konverguje kvadraticky. 11
Obr´ azek 6: Newtonova metoda pro ˇreˇsen´ı neline´arn´ı rovnice.
2.3
Newtonova metoda
Metoda bisekce nepouˇz´ıv´ a jin´e vlastnosti funkˇcn´ıch hodnot neˇz jejich znam´enka, coˇz vede k tomu, ˇze konverguje vˇzdy, ale pomalu. Pokud se vyuˇzij´ı tak´e velikosti funkˇcn´ıch hodnot, m˚ uˇzeme odvodit rychleji konverguj´ıc´ı metody, kter´e n´am v kaˇzd´e iteraci budou d´avat pˇresnˇejˇs´ı aproximaci koˇrene ˇreˇsen´e rovnice. V prvn´ı ˇradˇe se zde vyuˇz´ıv´a aproximace funkce f vystupuj´ıc´ı v rovnici pomoc´ı prvn´ıch dvou ˇclen˚ u jej´ıho Taylorova rozvoje, tedy f (x + h) ≈ f (x) + f 0 (x)h, coˇz je line´arn´ı funkce h, kter´ a aproximuje f v okol´ı bodu x. Nahrad´ıme tud´ıˇz neline´arn´ı funkci f touto line´ arn´ı funkc´ı, jej´ıˇz nulov´ y bod v h se snadno vypoˇc´ıt´a, je to h = −f (x)/f 0 (x), pokud 0 ovˇsem f (x) 6= 0. Je jasn´e, ˇze koˇreny obou tˇechto funkc´ı nejsou obecnˇe identick´e, takˇze popsan´ y postup mus´ıme iteraˇcnˇe opakovat. To vede k iteraˇcn´ı metodˇe, kter´e se ˇr´ık´a Newtonova metoda (nebo tak´e Newtonova-Raphsonova), jej´ıˇz algoritmus uv´ad´ıme jako Algoritmus 2. Algoritmus 2 Newtonova metoda x0 = poˇc´ ateˇcn´ı aproximace for k = 0, 1, 2, . . . xk+1 = xk − f (xk )/g 0 (xk ) end Ne obr´azku 6 ukazujeme, ˇze Newtonova metoda se d´a interpretovat jako aproximace funkce f pobl´ıˇz xk teˇcnou ke grafu funkce vedenou v bodˇe (xk , f (xk )). Jako dalˇs´ı aproximaci ˇreˇsen´ı pak bereme nulov´ y bod t´eto line´ arn´ı teˇcn´e funkce a proces postupnˇe opakujeme. Nˇekdy se Newtonovˇe metodˇe proto tak´e ˇr´ık´ a metoda teˇcen. Pˇr´ıklad 11 (Newtonova metoda). Newtonovu metodu pˇredvedeme opˇet na hled´an´ı koˇrene rovnice f (x) = x2 − 4 sin x = 0. Derivace t´eto funkce je f 0 (x) = 2x − 4 cos x, takˇze iteraˇcn´ı sch´ema je d´ ano vzorcem xk+1 = xk −
x2k − 4 sin xk . 2xk − 4 cos xk 12
Jako poˇc´ateˇcn´ı pˇribl´ıˇzen´ı zvol´ıme x0 = 3 a postupnˇe obdrˇz´ıme posloupnost iterac´ı, kter´a je uvedena d´ale. Pˇritom hk = −f (xk )/(xk ) oznaˇcuje zmˇenu xk v kaˇzd´e iteraci. Iteraˇcn´ı proces m˚ uˇzeme ukonˇcit, kdyˇz bude |hk |/|xk | nebo |f (xk )|, nebo oboj´ı, menˇs´ı neˇz n´ami pˇredepsan´ a tolerance. k
xk
f (xk )
f 0 (xk )
hk
0 1 2 3 4
3.000000 2.153058 1.954039 1.933972 1.933754
8.435520 1.294772 0.108438 0.001152 0.000000
9.959970 6.505771 5.403795 5.288919 5.287670
-0.846942 -0.199019 -0.020067 -0.000218 0.000000
Na Newtonovu metodu se m˚ uˇzeme tak´e d´ıvat jako na speci´aln´ı zp˚ usob pˇrevodu neline´arn´ı rovnice f (x) = 0 na u ´lohu o pevn´em bodˇe pro jistou funkci g, tedy x = g(x), kde za funkci g vol´ıme g(x) = x − f (x)/f 0 (x). a pevn´ y bod hled´ ame metodou postupn´ ych aproximac´ı. Abychom vyˇsetˇrili konvergenci metody, potˇrebujeme tedy nejprve zn´ at derivaci funkce g, coˇz je po u ´pravˇe g 0 (x) = f (x)f 00 (x)/(f 0 (x))2 (pokud f 0 (x) 6= 0). Je-li tedy x∗ jednoduch´ y koˇren, tj. f (x∗ ) = 0 a f 0 (x∗ ) 6= 0, pak g 0 (x∗ ) = 0. Newtonova metoda m´ a tedy pro jednoduch´e koˇreny asymptoticky kvadratickou rychlost konvergence, tedy r = 2. Kvadratick´a rychlost konvergence Newtonovy metody znamen´a, ˇze asymptoticky (v bl´ızkosti koˇrene) se chyba metody po kaˇzd´e iteraci umocn´ı na druhou. Jinak tak´e m˚ uˇzeme ˇr´ıci, ˇze se poˇcet pˇresn´ ych (spr´ avn´ ych) ˇc´ıslic pˇribliˇzn´eho ˇreˇsen´ı po kaˇzd´e iteraci zdvojn´asob´ı. Naproti tomu pro n´asobn´e koˇreny je Newtonova metoda pouze line´arnˇe (lok´alnˇe) konvergentn´ı s konstantou C = 1−(1/m), kde m je n´ asobnost poˇc´ıtan´eho koˇrene. Opˇet ale mus´ıme zd˚ uraznit, ˇze tyto u ´vahy o konvergenci plat´ı pouze lok´ alnˇe v nˇejak´em vˇetˇs´ım nebo menˇs´ım okol´ı hledan´eho koˇrene a ˇze Newtonova metoda, kter´ a nen´ı odstartov´ana dostateˇcnˇe bl´ızko ke koˇreni, nemus´ı konvergovat v˚ ubec. Jednoduch´ y pˇr´ıklad je situace, kdy nˇekdy bˇehem iterac´ı bude f 0 (xk ) relativnˇe mal´e (graf funkce f bude m´ıt v bodˇe xk t´emˇeˇr vodorovnou teˇcnu) a v d˚ usledku toho bude n´asleduj´ıc´ı iterace m´ıt tendenci leˇzet nˇekde daleko od posledn´ıho pˇribl´ıˇzen´ı. Pˇr´ıklad 12 (Newtonova metoda pro n´ asobn´ y koˇren). N´asleduj´ıc´ı dva pˇr´ıklady ukazuj´ı oba typy v´ yˇse popsan´eho chov´ an´ı Newtonovy metody. Prvn´ı z nich ukazuje kvadratickou konvergenci k jednoduch´emu koˇrenu, druh´ y line´ arn´ı konvergenci k n´asobn´emu koˇrenu. N´asobnost koˇrene ve druh´em z uveden´ ych pˇr´ıklad˚ u je 2, takˇze C = 1/2.
13
k 0 1 2 3 4 5
2.4
f (x) = x2 − 1
f (x) = x2 − 2x + 1
xk 2.0 1.25 1.025 1.0003 1.00000005 1.0
xk 2.0 1.5 1.25 1.125 1.0625 1.03125
Metoda seˇ cen
Jistou nev´ yhodou Newtonovy metody je, ˇze za jej´ı kvadratickou konvergenci plat´ıme t´ım, ˇze v kaˇzd´em iteraˇcn´ım kroku mus´ıme kromˇe funkˇcn´ı hodnoty poˇc´ıtat tak´e hodnotu derivace. V´ ypoˇcet hodnot derivace pˇritom m˚ uˇze b´ yt nepohodln´ y nebo ˇcasovˇe n´aroˇcn´ y, takˇze bychom mohli uvaˇzovat o tom, ˇze hodnoty derivac´ı budeme nahrazovat diferenˇcn´ımi pod´ıly vypl´ yvaj´ıc´ımi z definice derivace funkce, tedy bychom mohli pro vhodn´e dostateˇcnˇe mal´e h kl´ast napˇr´ıklad f 0 (x) ≈
f (x + h) − f (x) . h
To by ovˇsem znamenalo poˇc´ıtat v kaˇzd´e iteraci jednu funkˇcn´ı hodnotu nav´ıc, a to jen proto, abychom z´ıskali pˇribliˇznou informaci o hodnotˇe derivace. Lepˇs´ı je zaloˇzit podobnou diferenˇcn´ı aproximaci derivace na funkˇcn´ıch hodnot´ach, kter´e jsme uˇz bˇehem iterac´ı stejnˇe vypoˇc´ıtali, a kl´ast f (xk ) − f (xk−1 ) f 0 (xk ) ≈ . xk − xk−1 Tento postup vede k metodˇe seˇcen, jej´ıˇz algoritmus uv´ad´ıme jako Algoritmus 3. Na obr´azku 7 vid´ıme, ˇze metoda seˇcen se d´ a interpretovat jako aproximov´an´ı funkce f pˇr´ımkou proch´azej´ıc´ı pˇredchoz´ımi dvˇema iteracemi, tedy seˇcnou, pˇriˇcemˇz za nov´e pˇribl´ıˇzen´ı bereme nulov´ y bod t´eto line´arn´ı funkce. Na rozd´ıl od Newtonovy metody zde ovˇsem potˇrebujeme dvˇe poˇc´ateˇcn´ı aproximace. Algoritmus 3 Metoda seˇcen x0 , x1 = poˇc´ ateˇcn´ı aproximace for k = 0, 1, 2, . . . xk+1 = xk − f (xk )(xk − xk−1 )/(f (xk ) − f (xk − +)) end Pˇr´ıklad 13 (Metoda seˇcen). Metodu seˇcen budeme ilustrovat opˇet na hled´an´ı koˇrene rovnice f (x) = x2 − 4 sin x = 0. Za potˇrebn´a dvˇe poˇc´ ateˇcn´ı pˇribl´ıˇzen´ı vezmeme x0 = 1 a x1 = 3, vypoˇc´ıt´ame pˇr´ısluˇsn´e funkˇcn´ı hodnoty a za dalˇs´ı pˇribliˇzn´e ˇreˇsen´ı vezmeme pr˚ useˇc´ık pˇr´ımky spojuj´ıc´ı tyto dvˇe funkˇcn´ı hodnoty s nulou. Cel´ y postup pak opakujeme, pˇriˇcemˇz pouˇzijeme toto novˇe z´ıskan´e pˇribl´ıˇzen´ı koˇrene a tu novˇejˇs´ı ze dvou pˇredch´ azej´ıch iterac´ı, takˇze v kaˇzd´em iteraˇcn´ım kroku potˇrebujeme vypoˇc´ıtat pouze jednu novou funkˇcn´ı hodnotu. Posloupnost proveden´ ych iterac´ı je uvedena v tabulce, kde hk oznaˇcuje zmˇenu xk v pˇr´ısluˇsn´e iteraci. 14
Obr´ azek 7: Metoda seˇcen pro ˇreˇsen´ı neline´arn´ı rovnice.
k
xk
f (xk )
hk
0 1 2 3 4 5 6 7 8
1.000000 3.000000 1.438070 1.724805 2.029833 1.922044 1.933174 1.933757 1.933754
-2.365884 8.435520 -1.896774 -0.977706 0.534305 -0.061523 -0.003064 0.000019 0.000000
-1.561930 0.286735 0.305029 -0.107789 0.011130 0.000583 -0.000004 0.000000
Protoˇze kaˇzd´e nov´e pˇribliˇzn´e ˇreˇsen´ı, kter´e d´av´a metoda seˇcen, z´avis´ı na dvou pˇredchoz´ıch iterac´ıch, je vyˇsetˇrov´ an´ı konvergence metody o nˇeco sloˇzitˇejˇs´ı a detaily jsme nuceni zde vypustit. Uv´ad´ıme alespoˇ n, ˇze se d´ a dok´ azat, ˇze chyby metody splˇ nuj´ı pro jistou kladnou konstantu c > 0 vztah |ek+1 | lim = c, k→∞ |ek | · |ek−1 | coˇz znamen´a, ˇze posloupnost iterac´ı metodou seˇcen lok´alnˇe konverguje a rychlost konvergence je superline´arn´ı. Pˇresnˇeji (viz [1], [2] a [3]) se d´a uk´azat, ˇze asymptotick´a rychlost konvergence metody seˇcen je1 √ 1+ 5 ≈ 1,618. r= 2 Stejnˇe jako u Newtonovy metody je i u metody seˇcen ke konvergenci nutno iterace odstartovat dostateˇcnˇe bl´ızko koˇrene. Porovn´ame-li metodu seˇcen s Newtonovou metodou, vid´ıme, ˇze metoda seˇcen m´a v´ yhodu v tom, ˇze v kaˇzd´e iteraci potˇrebuje vypoˇc´ıtat pouze jednu novou funkˇcn´ı hodnotu. Za nev´ yhodu bychom mohli povaˇzovat to, ˇze vyˇzaduje dvˇe startovac´ı hodnoty a ˇze v˚ uˇci Newtonovˇe metodˇe konverguje pomaleji, i kdyˇz st´ ale superline´arnˇe. Menˇs´ı pracnost proveden´ı jedn´e iterace vyv´aˇz´ı u metody seˇcen zpravidla to, ˇze k dosaˇzen´ı koneˇcn´eho v´ ysledku mus´ıme prov´est vˇetˇs´ı poˇcet iterac´ı. D´a se tedy ˇr´ıci, ˇze nalezen´ı pˇribliˇzn´e hodnoty ˇreˇsen´ı neline´arn´ı rovnice metodou seˇcen je ˇcasto m´enˇe pracn´e neˇz pouˇzit´ı Newtonovy metody. 1
Pokud v´ am to ˇc´ıslo pˇripad´ a povˇedom´e, pˇripom´ın´ am, ˇze je to hodnota zlat´eho ˇrezu.
15
3 3.1
Dodatky Bezpeˇ cn´ e metody
Rychle konverguj´ıc´ı metody pro numerick´e ˇreˇsen´ı neline´arn´ıch rovnic jako jsou napˇr´ıklad Newtonova metoda ˇci metoda seˇcen (dalˇs´ı takov´e metody lze naj´ıt v literatuˇre [1, 2, 3]) nejsou bezpeˇcn´e v tom smyslu, ˇze pokud nejsou odstartov´any dostateˇcnˇe bl´ızko koˇrene, nemus´ı konvergovat. Bezpeˇcnou metodou v tomto smyslu je metoda p˚ ulen´ı intervalu, kter´a je ale pomal´ a a tedy n´akladn´ a. Jakou metodu tedy volit? ˇ sen´ım tohoto dilematu jsou hybridn´ı metody, kter´e jsou zahrnuty ve vˇetˇsinˇe modern´ıho maReˇ tematick´eho softwaru a kter´e v sobˇe kombinuj´ı vlastnosti obou v´ yˇse popsan´ ych typ˚ u metod. Jejich algoritmy jsou ale ovˇsem sloˇzitˇejˇs´ı. Tyto metody mohou napˇr´ıklad pracovat s rychle konvergentn´ı metodou a pˇritom doc´ılit toho, ˇze iterace z˚ ust´avaj´ı uvnitˇr poˇc´ateˇcn´ı uz´avˇery koˇrene. Pokud n´asleduj´ıc´ı aproximace ˇreˇsen´ı rychl´ ym algoritmem padne mimo interval uz´avˇery, vr´at´ıme se a provedeme jednu iteraci bezpeˇcnou metodou, napˇr´ıklad bisekc´ı. Pak se m˚ uˇze zkusit opˇet pouˇzit´a rychl´ a metoda, tentokr´ at ovˇsem uˇz na menˇs´ım intervalu a s vˇetˇs´ı nadˇej´ı na u ´spˇech. ke konci v´ ypoˇctu uˇz by mˇely iterace bˇeˇzet tou rychlou metodou. Uveden´ y postup je jen zˇr´ıdka horˇs´ı neˇz pouˇzit´ a pomal´ a metoda, zpravidla je mnohem rychlejˇs´ı. Popul´arn´ı implementace v´ yˇse popsan´eho hybridn´ıho postupu dnes poch´az´ı od Brenta (v literatuˇre tak´e tedy Brentova metoda) a kombinuje v sobˇe bezpeˇc´ı bisekce s rychlejˇs´ı konvergenc´ı tzv. inverzn´ı kvadratick´ e interpolace (v´ıce k tomu viz [3]). D´ıky tomu, ˇze se zde vyh´ yb´ame Newtonovˇe metodˇe, nejsou k v´ ypoˇctu zapotˇreb´ı hodnoty derivace. Soudob´ y kvalitn´ı software mus´ı pˇri implementaci metody vz´ıt v u ´vahu tak´e to, ˇze se jej´ı algoritmus realizuje v poˇc´ıtaˇcov´e aritmetice, tedy napˇr. ohl´ıdat moˇzn´ a pˇrekroˇcen´ı rozsahu poˇc´ıtaˇce nebo nepˇrimˇeˇrenˇe pˇr´ısn´e poˇzadavky na pˇresnost v´ ysledku. Dobrou implementaci v´ yˇse popsan´eho postupu pˇredstavuje napˇr´ıklad funkce fzero v Matlabu. Poznamen´av´ ame jeˇstˇe, ˇze jakousi kombinac´ı metody bisekce a metody seˇcen je metoda regula falsi (z lat., doslova pravidlo falˇse). Kaˇzd´ y jej´ı krok zaˇc´ın´a t´ım, ˇze body xk a xk−+ tvoˇr´ı uz´avˇeru hledan´eho koˇrene, ale m´ısto aby se v kaˇzd´em kroku interval uz´avˇery p˚ ulil, vypoˇc´ıt´ a se nejprve xk+1 pomoc´ı vzorce metody seˇcen. Pr˚ ubˇeh funkce se tedy na dan´em intervalu opˇet nahrad´ı seˇcnou. Pak se z takto z´ıskan´ ych tˇr´ı bod˚ u zachovaj´ı ty dva, v nichˇz m´a funkce f opaˇcn´ a znam´enka, a postup se opakuje. Metoda regula falsi je dalˇs´ı vˇzdy konvergentn´ı metodou, mus´ıme ji ovˇsem odstartovat z uz´ avˇery koˇrene. Jej´ı konvergence je pouze line´arn´ı a m˚ uˇze, ale nemus´ı b´ yt rychlejˇs´ı neˇz metoda p˚ ulen´ı. Lze tak´e uk´azat, ˇze v nˇekter´ ych pˇr´ıpadech m˚ uˇze jeden z krajn´ıch bod˚ u uz´avˇery z˚ ust´ avat bˇehem iterac´ı trvale beze zmˇeny a aˇckoli druh´ y bod konverguje ke koˇrenu rovnice, uz´avˇera se nem˚ uˇze zmenˇsit pod jistou mez.
3.2
Numerick´ y v´ ypoˇ cet koˇ ren˚ u polynomu
Aˇz dosud jsme se zab´ yvali metodami pro nalezen´ı jednoho nulov´eho bodu obecn´e re´aln´e funkce jedn´e re´aln´e promˇenn´e. Pokud je uvaˇzovan´a funkce polynom p(x) stupnˇe n, pak potˇrebujeme ˇcasto naj´ıt vˇsechny jeho nulov´e body, z nichˇz nˇekter´e mohou b´ yt komplexn´ı, i kdyˇz polynom s´am m´a re´aln´e koeficienty. O koˇrenech polynom˚ u n´am algebraick´a teorie ˇr´ık´a podrobnˇejˇs´ı informace neˇz zn´ame o nulov´ ych bodech obecn´ ych funkc´ı. Pˇredevˇs´ım je zde tzv. z´ akladn´ı vˇeta algebry, podle n´ıˇz kaˇzd´ y polynom stupnˇe n m´ a v komplexn´ı rovinˇe pr´avˇe n nulov´ ych bod˚ u (koˇren˚ u), pokud kaˇzd´ y z nich poˇc´ıt´ ame tolikr´ at, kolik ˇcin´ı jeho n´asobnost. D´ale se d´a uk´azat, ˇze pokud m´a re´aln´ y polynom komplexn´ı koˇreny, vyskytuj´ı se tyto koˇreny vˇzdy ve dvojic´ıch komplexnˇe 16
sdruˇzen´ ych ˇc´ısel, tedy jako x ± ıy. Pro hled´an´ı koˇren˚ u polynom˚ u nen´ı nezbytn´e pouˇz´ıvat komplexn´ı aritmetiku, leckdy lze poˇc´ıtat jejich re´aln´e a imagin´ arn´ı ˇc´ asti x a y oddˇelenˇe. Pro v´ ypoˇcet koˇren˚ u polynom˚ u existuje ˇrada moˇznost´ı: • Pouˇzijeme nˇekterou z popsan´ ych obecn´ ych metod(napˇr. Newtonovu metodu) a nalezneme jeden koˇren x1 . Pak d´ ale pracujeme s redukovan´ ym polynomem p(x)/(x−x1 ), jehoˇz stupeˇ n je o jedniˇcku niˇzˇs´ı. Postup opakujeme tak dlouho, dokud nestanov´ıme vˇsechny koˇreny. Metoda se komplikuje, pokud naraz´ıme na komplexn´ı koˇren. • K dan´emu polynomu sestav´ıme jeho doprovodnou matici, coˇz je speci´aln´ı matice maj´ıc´ı vlastn´ı ˇc´ısla shodn´ a s koˇreny polynomu. Pak nˇejakou vhodnou numerickou metodou algebry stanov´ıme jako koˇreny dan´eho polynomu vlastn´ı ˇc´ısla t´eto matice. Tento postup, kter´ y je pouˇzit ve funkci roots v Matlabu, je spolehliv´ y, ale nen´ı tak efektivn´ı jako pouˇzit´ı numerick´ ych metod odvozen´ ych speci´alnˇe pro v´ ypoˇcet koˇren˚ u polynomu. • Pouˇzijeme nˇekterou ze speci´ aln´ıch metod pro v´ ypoˇcet nulov´ ych mod˚ u polynom˚ u. Najdou se mezi nimi jak bezpeˇcn´e metody, kter´e izoluj´ı koˇreny napˇr´ıklad ve sjednoceni disk˚ uv komplexn´ı rovinˇe (ty jsou ovˇsem podobnˇe jako bisekce pouze line´arnˇe konvergentn´ı), tak rychle konverguj´ıc´ı metody (i rychlejˇs´ı neˇz je Newtonova metoda). O tˇechto speci´aln´ıch metod´ ach se lze pouˇcit napˇr´ıklad v [1, 2].
3.3
Numerick´ eˇ reˇ sen´ı soustav neline´ arn´ıch rovnic
ˇ sen´ı soustav neline´ Reˇ arn´ıch rovnic je obt´ıˇznˇejˇs´ı, neˇz je tomu u jedn´e rovnice, a to z ˇrady d˚ uvod˚ u: • Chov´an´ı soustavy m˚ uˇze b´ yt mnohem rozmanitˇejˇs´ı neˇz chov´an´ı jedn´e rovnice (a jejich koˇren˚ u). Teoretick´ a anal´ yza existence a poˇctu ˇreˇsen´ı je tak mnohem sloˇzitˇejˇs´ı. • Konvenˇcn´ı metody pouˇz´ıvan´e pro jednu rovnici se leckdy daj´ı v´ıcem´enˇe pˇr´ımoˇcaˇre zobecnit i pro soustavy, ale u soustav nen´ı jednoduch´ y zp˚ usob, jak zobecnit pojem uz´avˇery ˇreˇsen´ı, takˇze zde nen´ı jednoduch´e sestrojit bezpeˇcn´e, glob´alnˇe konverguj´ıc´ı metody. Urˇcit´e moˇznosti zde ale existuj´ı, nicm´enˇe se vymykaj´ı moˇznostem tohoto textu a nenajdou se ani v bˇeˇzn´ ych uˇcebnic´ıch. Nicm´enˇe v Matlabu je pro ˇreˇsen´ı soustav neline´arn´ıch rovnic k dispozici vcelku spolehliv´ a funkce fsolve. • Pracnost numerick´eho ˇreˇsen´ı soustav neline´arn´ıch rovnic roste neline´arnˇe s poˇctem nezn´am´ ych. Tak napˇr´ıklad jeden iteraˇcn´ı krok Newtonovy metody pro soustavu o n nezn´am´ ych znamen´a obecnˇe v´ ypoˇcet n2 hodnot derivac´ı a jedno ˇreˇsen´ı souatavy n line´arn´ıch rovnic o n nezn´am´ ych, coˇz je samo o sobˇe obecnˇe ˇr´adovˇe n3 aritmetick´ ych operac´ı. Tak´e organizaˇcn´ı struktura algoritm˚ u je mnohem sloˇzitˇejˇs´ı. Jak jsme uvedli jiˇz na zaˇc´ atku tohoto textu, studium problematiky soustav se zde vymyk´a naˇsim ´ ˇcasov´ ym moˇznostem. Uvodn´ ı informace m˚ uˇze z´ajemce naj´ıt v doporuˇcen´e literatuˇre [1], [2], [3].
Reference ˇ Numerick´e metody algebry. Praha: SNTL, 1982. [1] M´IKA, Stanislav. MVST. 17
ˇ 2002. [2] M´IKA, Stanislav a Marek BRANDNER. Numerick´e metody I. Plzeˇ n: FAV ZCU, [3] HEATH, M. T. Scientific Computing: An Introductory Survey. 2nd Edition. New York: McGraw-Hill, 2002, 563 s.
18
