Numerické metody — nalezení kořene funkce metodou bisekce

Proč numerické metody?

Někdy se neskutečně hodí znát řešení[1]kořeny nějaké rovnice n-tého řádu, ale není třeba znát naprosto přesné analytické řešení, ale řekněme nějaké řešení s přesností na \(3\) platné číslice. Proč bychom vůbec ale takové řešení chtěli?

Možné důvody, proč něco takového chceme

Jedním z důvodů takového nepřesného řešení může být např. to, že chceme přibližný výsledek zobrazit na nějakém displeji s určitým počtem číslic, případně pro nějaké měření či výstupy je nám poměrně k ničemu vědět, že je trubka dlouhá \(\sqrt{3}\ \mathrm{m}\), ale prostě stačí vědět, že je zhruba \(1,7\ \mathrm{m}\).

Důvodem takové nepřesnosti může být občas prostě jen to, že to takto stačí, anebo také, že ne vždy je možnost nějaké analytické řešení takto snadno najít.

Práce s chybou

Při práci s numerickými metodami téměř vždy[2]Výjimkou budiž řešení, kde se do přesného řešení prostě náhodou trefíme. budeme pracovat s nějakou chybou. Tato chyba je způsobena samotným principem toho, jak na daná řešení přicházíme — využíváme totiž nějakého (pokud možno) sofistikovaného odhadu. Ukažme si to třeba na příkladu, mějme např. takovouto rovnici:

$$x^3 – 3 = 0$$

Je asi zcela evidentní, že pro řešení v \(x\in \mathbb{R}\) hodnota bude \(x = \sqrt[3]{3}\), což je velmi přibližně \(1,44\)[3]1,4422495703074083823216383107801095883918692534993505775464161945…. Nicméně víme, že použijeme-li hodnotu \(1,44\), pracujeme s chybou[4]resp. absolutní chybou \(\mathrm{e}=\mid{1,44-\sqrt[3]{3}}\mid\).  Otázkou tedy jest, zda-li nám takováto přesnost stačí. Pro uříznutí trubky zubatou pilou zcela jistě ano, pro nastavení dorazu CNC frézy pro výbrus válce už bych si tak jistý nebyl.

Metoda bisekce

Metoda bisekce využívá faktu, že alespoň trochu známe průběh funkce, jejíž kořeny (tedy body, pro které platí \({f}(x)=0\)) hledáme. Tato funkce musí být v intervalu zkoumání spojitá. Podíváme-li se na průběh výše uvedené funkce \(f(x): y = x^3-3\), dostaneme zhruba:

Graf funkce: \(y=x^3-3\)
Graf funkce: \(y=x^3-3\)

Vidíme, že funkce protne osu \(x\) opravdu někde v intervalu \(x \in (1; 1,5)\) Čeho si u takové funkce můžeme všimnout? Vezmeme-li nějaký bod \(x_0\), který bude naším opravdovým přesným kořenem (tedy \(\sqrt[3]{3}\),) poté určitě bude platit, že součin \(f(x_0-x_e) \cdot f(x_0 + x_e) \leq 0\) pro \( x_e \in \mathbb{R}, x_e \geq 0 \). Toto vidíme krásně z grafu, proč platí. Vezmu-li nějaké dvě záporné hodnoty a vynásobím je, dostanu kladné číslo. A protože vidím, že funkce zde osu \(x\) protíná, tedy nejdříve je v záporné části a pak v kladné, vezmu-li dva body kolem tohoto bodu průniku (kořene), je vidět, že nutně musím vybrat jeden kladný a jeden záporný bod. A tedy součin dvou funkčních hodnot takovéto funkce kolem kořene rovnice bude vždy záporný.

Metoda bisekce tohoto faktu využívá, tedy vytvoříme odhad \(x_s\) na hodnotu našeho kořene \(x_0\) a zjišťujeme (po dosazení tutoho odhadu), jak moc se lišíme oproti nulové hodnotě. Označím-li spodní (dolní) odhad \(x_d\) a horní odhad jako \(x_h\) a takové, že \( (x_d < x_0) \wedge (x_h > x_0)\), resp. tedy \(x_d < x_0 < x_h\), pak se zaručeně musíme body z takového intervalu trefit do jedné ze tří možností:

  1. Trefíme náhodně přesné řešení, tedy např. kdy \(f(x_d)=0\) a tedy \(x_d = x_0\). Poté algoritmus výpočtu můžeme úspěšně ukončit. Samozřejmě toto přesné řešení nemůžete trefit, pokud kořen není přesně vyjádřitelný, tedy s ukončeným desetinným rozvojem.
  2. Nacházíme se vedle přesného řešení vlevo, tedy \(x_s < x_0\) a součin \(f(x_d)\cdot f(x_s) > 0\) — jak jsem psal výše, budeme-li násobit dvě záporná čísla mezi sebou, výsledek tohoto součinu bude jistě kladný.
  3. Nacházíme se vedle přesného řešení vpravo, tedy \(x_s > x_0\) a součin \(f(x_d)\cdot f(x_s) < 0\), tedy podle podobné logiky jako v druhém bodě, zde budu-li násobit jedno záporné číslo s kladným, výsledek součinu bude určitě záporný.

Jak tento postup algoritmizovat

Metoda bisekce je neskutečně jednoduchá, má však svoje zásadní úskalí — vždy je její přesnost závislá na několika faktorech: jak dobře uděláme odhad chování průběhu funkce kolem kořene a jak moc dobře odhadneme množství kořenů. Bisekce totiž dokáže pracovat s právě jedním kořenem naráz a pro různé kořeny je potřeba počítat zvlášť. Ukažme si však nyní na tomto případu rovnice \(x^3-3=0\), jak metodu bisekce použít, případně naprogramovat.

Odhadnutí průběhu funkce

Zde tak nějak vidíme, že funkce bude \(y=x^3-3\) bude rostoucí a z grafu funkce vidíme, že rovnice \(0 = x^3-3\) bude mít zcela jistě pouze jediný kořen v \(\mathbb{R}\)[5]Ve skutečnosti má kořeny 3 v \(\mathbb{C}\), akorát jeden kořen má nulovou imaginární složku a dva zbylé nikoliv. Protože však funkci řešíme pouze v \(\mathbb{R}\), bude nás zajímat … Continue reading. Nyní je tedy čas na odhad tohoto kořene. U takto jednoduché funkce, kdy víme, jak vypadá průběh funkce \(y=x^3\), kterou posuneme o \(3\) jednotky dolů máme zhruba představu, že kořen se bude nacházet určitě někde mezi řekněme jedničkou a trojkou. Ne vždy však takový fundovaný odhad budeme schopni od oka udělat, funkce může mít například více kořenů a její průběh může být prostě velmi různorodý.

V takovém případě pomůže jednoduchá tabulka, prostě zkusíme pár hodnot a uvidíme, kde jsme přelezli nulu. Samozřejmě, že takový postup je citlivý na nastavení vhodného kroku, protože pokud krok bude příliš velký, můžeme nulu přelézt třeba i několikrát a vůbec si nevšimneme. Pokud bude krok příliš malý, budeme počítat dlouho (případně hodně dlouho.) V případě počítání takto jednoducých funkcí je to prakticky jedno, ale pokud budeme počítat např. mesh (tedy nějakou síť) v rozlišení \(100×100\) uzlů, budeme počítat paralelně \(10\,000\) rovnic a pokud každou zesložitíme třeba desetkrát, už to prostě na délce výpočtu bude znát.

Horní a dolní mez

Vytvořme však nyní pro naši funkci \(y=x^3-3\)  malou tabulku s několika dosazenými hodnotami:

x-2-10123
y-11-4-3-2524

V tabulce vidíme, že osu \(x\) funkčními hodnotami přelezeme někde v intervalu \(x\in(1; 2)\), tedy tyto hodnoty v klidu můžeme použít jako výchozí krajní odhady pro náš výpočet.

Výpočet \(x_s\) a jeho užití

Prohlásíme tedy, že \(x_d=1, x_h=2\) a toto bude náš první krok. Nyní je na čase odhadnout nějakým způsobem bod uprostřed, který vybereme. Máme vícero možností, jak na to jít, ale nejpřímočařejší metodou je prosté rozpůlení intervalu, tedy najdeme takové číslo, které se nachází mezi horní a dolní mezí přesně uprostřed — jinými slovy najdeme klasický aritmetický průměr, tedy \(x_s = \frac{x_h + x_d}{2} = \frac{1 + 2}{2} = \frac{3}{2} = 1,5\).

Nyní musíme rozhodnout, jak dobře jsme se trefili, čili zhodnotit, na jaké straně odhadu oproti přesné hodnotě jsme. Půjdeme na to čistě logicky — je jasné, že pokud vynásobíme mezi sebou dvě záporné či dvě kladné hodnoty, dostaneme kladnou, pokud kladnou a zápornou, dostaneme zápornou hodnotu. Této znalosti snadno využijeme, tedy víme, že pokud vynásobíme hodnotu funkce ve spodním limitu \(x_d\), která bude pro rostoucí funkci určitě záporná, pak víme, že pokud jsme se trefili pod přesnou hodnotu, tedy \(x_s<x_0\), bude i funkční hodnota \(f(x_s)\) pro náš bod záporná a součin těchto dvou čísel bude kladný a samozřejmě naopak pro trefu nad přesnou hodnotu, tedy \(x_s>x_0\). Případně v takto jednoduchém případě můžeme rovnou hledat pouze kladnost či zápornost \(f(x_s)\).

Zpřesňování výsledku

A nyní rozhodovací možnost: pokud jsme se s naší střelbou trefili tak, že \(x_s<x_d\), poté víme určitě, že řešení se nachází někde pro číslo větší než \(x_s\). Proto můžeme bezpečně prohlásit náš nový spodní limit \(x_d\) právě číslo \(x_s\). Celé tyto kroky zopakujeme, budeme neustále hledat průměrné hodnoty mezi \(x_d\) a \(x_h\) a horní či dolní hranici posouvat tak, abychom postupně zužovali místo, kde se daný kořen nachází. Nasťiňme si nyní jednotlivé kroky, jak bude výpočet probíhat:

  1. Dolní mez nastavíme na \(x_d = 1\), horní mez na \(x_h)=2\). Průměrná hodnota mezi nimi je \(x_s=1,5\).  Dosadíme tuto průměrnou hodnotu do samotné funkce a zjistíme, že \(f(x_s) = f(1,5) = 0,375\), tedy kladné číslo. Zcela jistě jsme se tak trefili někam nad hodnotu \(x_0\), kterou hledáme.
  2. Dolní mez tedy necháme na \(x_d=1\), ale horní posuneme na nám známou hodnotu \(x_h = x_s = 1,5\). Najdeme střední hodnotu \(x_s = \frac{1+1,5}{2} = \frac{2,5}{2} = 1,25\). Tuto hodnotu dosadíme zpět do původní funkce a zjistíme, jak se v tomto bodě chová: \(f(x_s) = f(1,25) \doteq -1\), tedy  výsledek je záporný a víme, že jsme svým odhadem překročili hodnotu \(x_0\). Nyní tedy máme jistotu, že výsledek se nachází v intervalu \(x_0 \in (1,25; 1,5)\).
  3. Nastavme tedy tuty hodnoty opět jako naše meze a zkoumejme dále, jak se bude chovat rozpůlený interval, tedy \(x_d=1,25, x_h=1,5, x_s = \frac{1,25}{1,5} = 1,375\) a tedy \(f(x_s) \doteq -0,4\) a tedy záporná hodnota.  Stále jsme tedy se svým odhadem \(x_s\) pod přesným řešením a můžeme opět zpřesnit dolní mez výpočtu. Nyní tedy jistě víme, že vysledek bude v intervalu \(x_0 \in (1,375; 1,5) \).

Tento postup je možné opakovat tak dlouho, dokud nás to baví, ale samozřejmě dříve či později musí nutně skončit. Pokud nám stačí takováto přesnost, můžeme klidně prohlásit, že výsledek je v tomto intervalu. Chceme-li však přesnější výpočet, musíme s bisekcí pokračovat.

Kdy zastavit výpočet?

Je asi jasné, že pro přesný výpočet by takovýto výpočet trval do nekonečna[6]Resp. do timeoutu či do vyčerpání počítatele., takže je třeba algoritmus nějakým způsobem omezit, aby výsledek byl konečný, byť nepřesný. Máme tyto zastavovací podmínky (musí platit alespoň jedna z nich):

  1. Trefíme se přesně, tedy \(f(x_d) = 0 \vee f(x_h)=0\). Poté máme přesné řešení a netřeba dále počítat (za řešení prohlásíme tu hodnotu, která dává funkční hodnotu nulovou.
  2. Výpočet trvá moc dlouho, tedy můžeme nastavit maximální počet kroků, které se mají dělat. Pokud se v potřebném počtu kroků nepovede spočítat výsledek na potřebnou přesnost (viz níže,) výpočet je ukončen s nepřesným výsledkem, resp. s chybou větší než požadovanou.
  3. Výsledek se nachází v rámci požadované chyby, tedy \(x_h-h_d < e\), kde \(e\) je chtěná přesnost (či nepřesnost, záleží na úhlu pohledu.)

Vyzkoušejte sami

Zdrojový kód k tomuto výpočtu zde:bisekce.zip
Případně můžete přejít přímo na stránku skriptu: Bisekce (skript).

Za konzultaci a proofreading děkuji Lukášovi a Mirouškovi 🙂

Poznámky pod čarou

Poznámky pod čarou
1 kořeny
2 Výjimkou budiž řešení, kde se do přesného řešení prostě náhodou trefíme.
3 1,4422495703074083823216383107801095883918692534993505775464161945…
4 resp. absolutní chybou
5 Ve skutečnosti má kořeny 3 v \(\mathbb{C}\), akorát jeden kořen má nulovou imaginární složku a dva zbylé nikoliv. Protože však funkci řešíme pouze v \(\mathbb{R}\), bude nás zajímat pouze kořen s nulovou imaginární složkou.
6 Resp. do timeoutu či do vyčerpání počítatele.

Napsat komentář

Vaše e-mailová adresa nebude zveřejněna. Vyžadované informace jsou označeny *