TI-Basic - Bisekce

Doufám, že dobrým úvodem do programování na kalkulačkách pro Vás bude realizace jednoho z nejjednodušších matematických algoritmů. Nečekejte žádnou efektivitu kódu či velkou rychlost výpočtu. Jde opravdu jen o ukázku realizace algoritmu, jelikož kalkulačky TI mají implementovány daleko sofistikovanější a rychlejší řešení tohoto problému.....

Teoretický úvod

Bisekce je jedna z nejzákladnějších, nejjednodušších, ale také nejnepřesnějších numerických metod pro zjišťování kořenů funkcí - jejich nulových bodů. Základní problém vyjádřený slovně a matematicky by zněl asi takhle:

Hledáme číslo x, pro které je funkční hodnota dané funkce rovna nule f(x)=0

Tato úloha se dá řešit analyticky i numericky. První způsob nabízí velkou výhodu přesného výsledku, ale v prostředí jednoduchých procesorů se dá realizovat těžko. Proto bylo vymyšleno mnoho způsobů řešení numerických a jedno z nich je právě řešení pomocí bisekce. Tato metoda je pro nás pouze instruktážní a bude sloužit pouze za účelem seznámení se s programováním v TI-Basicu. Kalkulačka nám samozřejmě nabízí mnoho způsobů, jak podobný problém vyřešit a svou přesností a rychlostí bisekci daleko předčí. I tak si ale myslím, že mnoho postupů uvedených zde se vám může hodit.

Pro pochopení základního principu nám pomůže graf zkoumané funkce. Pro jednoduchost zvolím kvadratickou funkci:

f(x)=x^2-2, chcete-li y=x^2-2.

Je to graf posunuté paraboly

Metoda bisekce vychází z jednoduchého principu:

Zvolíme si dva body a, b na ose x tak, aby a bylo menší než b. Samozřejmě musí oba body ležet v definičním oboru naší funkce, tj. funkce musí mít v obou bodech smysl. Vypočítáme funkční hodnoty v obou bodech f(a), f(b)

f(a) = a^2-2

f(b) = b^2-2

Tyto hodnoty nám pomůžou najít řešení naší úlohy. Platí totiž jedna matematická věta, která zní:

Jestliže existuje funkce f(x), která je spojitá na intervalu <a, b> a platí nerovnost f(a)*f(b) < 0, potom existuje bod c z intervalu (a, b), pro který platí, že f(c)=0.

Slovy prostého člověka: hledáme takový interval <a, b>, v jehož jednom krajním bodě a je funkční hodnota naší funkce kladná a ve druhém bodě b záporná (nebo naopak), tedy součin a*b bude záporný. Aby funkce přešla z kladných hodnot do záporných (nebo naopak), musí projít osou x a to je právě bod, který je nulovým bodem funkce - funkční hodnota v tomto bodě je rovna 0.

První zádrhel nastane v okamžiku, kdy se zeptáme, kde vlastně máme body a a b hledat, protože může nastat případ, kdy má některá funkce více, nebo dokonce nekonečně mnoho, nulových bodů (kořenů). Obecně neexistuje algoritmus, který by je numerickou formou našel. Musíme použít náš matematický talent a body odhadnout. Z podmínky f(a)*f(b) < 0 můžeme svoji volbu ověřit, protože se už jedná o matematický a booleovský úkon. Jestliže je náš první odhad správný, můžeme se přibližovat k nulovému bodu, protože máme interval, kde se nachází. Název bisekce znamená "dvě části". Budeme totiž dělit náš interval na dvě části a pro jednoduchost na dvě části stejné. Střed intervalu < a,b > nazveme s1. Je zřejmé, že náš nulový bod se nachází v některém ze dvou intervalů (a,s1) nebo (s1,b). Zjistíme to opět pomocí podmínky existence nulového bodu na intervalu přenesenou na naše nové intervaly: f(a)*f(s1) < 0 a f(s1)*f(b) < 0. Interval, který danou nerovnost splňuje, je náš nový favorit, protože v něm se nulový bod nachází. Hledáme dál: dejme tomu, že podmínku splnil první interval (a,s1). rozdělíme ho opět na dva stejné intervaly (a,s2) a (s2,s1) a otestujeme, ve kterém že se to ten nulový bod nachází. Stejný postup budeme aplikovat tak dlouho, až se trefíme na nulový bod. A jak to poznáme? Pro nulový bod platí vztah f(c)=0. Až se tedy bude součin f(a)*f(b) = 0, bude jeden z bodů a,b nulovým bodem funkce. Toto se může stát i docela rychle - záleží na volbě intervalu, ale taky velmi pomalu - to už záleží na charakteru nulového bodu. V případě, že nulový bod je číslo s neukončeným desetinným rozvojem, výsledku se nedočkáme, resp. narazíme na možnosti výpočetní techniky - citlivost (přesnost nebo rozsah). Proto si přesnost výpočtu hned na začátku zvolíme. Pokud vyjádříme chybu jako

d = (b-a)/2^k

kde a a b jsou hranice intervalu a k je konečný počet přiblížení (iterací) ke kořenu, můžeme porovnat výslednou chybu s maximální požadovanou a pokud je výsledná chyba menší, budeme považovat jeden z bodů konečného intervalu za dostatečně přesný nulový bod. Který bod onoho intervalu to bude, je v podstatě jedno, protože daný vztah pro výpočet chyby je vlastně délkou tohoto intervalu. Výsledek ve tvaru bod ± chyba tedy vyhovuje našim požadavkům.

Pojďme si tedy ukázat dělení intervalu na příkladu f(x)=x^2-2 s požadovanou maximální chybou d=0.1

Zvolili jsme interval (a,b), kde a=1 a b=2. Jejich funkční hodnoty jsou f(1) = -1 a f(2) = 2. Je zřejmé, že nerovnost f(1)*f(2) < 0 je splněna, budeme tedy tento interval dělit na polovinu. Střed spočítáme jako aritmetický průměr hodnot a a b:

s1 = (a+b)/2 = (1+2)/2 = 1.5

Vypočítáme funkční hodnotu v bodě s1 a ověříme nerovnosti v obou nových intervalech (a,s1) a (s1,b), tedy (1,1.5) a (1.5,2). Nerovnost f(1.5)*f(2) < 0 neplatí, budeme tedy hledat dál v prvním intervalu. Pro ověření, zda už jsme náhodou nesplnili přesnost, dosadíme do vztahu pro chybu d = (s1-a)/2^1 = 0.25 a porovnáme s požadovanou přesností (d=0.1). Přesnost ještě není dostatečná, rozdělíme tedy tento interval na dva intervaly (a,s2) a (s2,s1). Otestujeme nerovnost a zjistíme, že nerovnost platí pro druhý interval (s2,s1). Opět otestujeme chybu a jak se můžeme přesvědčit, přesnost se nám dostala pod stanovenou hranici 0.1:

d = (s2-s1)/2^2 = 0.0625.

V případě, že bychom požadovali přesnost větší, hledáme stejným způsobem dál.

Základní algoritmus

1) požadujeme vstupní data: funkci f(x), počáteční interval (a,b), kde a < b a požadovanou chybu výpočtu (přesnost)

2)Otestujeme f(a)*f(b)

f(a)*f(b)=0 » určíme, která ze dvou hodnot a,b je nulovým bodem funkce, chyba je rovna 0, ukončíme program.

f(a)*f(b) < 0 » vypočítáme chybu a porovnáme s požadovanou chybou. Jestliže vyhovuje, ukončíme program

» rozdělíme interval (a,b) na dva podintervaly (a,s) a (s,b) a každý otestujeme f(a)*f(s) < 0 a f(s)*f(b) < 0

» interval, který podmínku splňuje označíme jako interval (a,b) a vrátíme se k bodu 2)

Algoritmus je již optimalizován, je tedy určený přímo k přepsání do programovacího jazyku. Optimalizace spočívá v nahrazení starého a již ověřeného intervalu (a,b) intervalem novým, který vyhověl nerovnosti. Odpadá tak nutnost ukládat každý nový střed intervalu do nové proměnné. Jednoduše nazveme střed intervalu novým áčkem nebo béčkem a k dalšímu ověření tohoto intervalu můžeme použít stejný algoritmus bez změny názvů proměnných.

Algoritmus programu bisekce()

Deklarace proměnných, nastavení režimu kalkulačky, vymazání I/O,

2. Vstup: Zadání funkce f(x)

3. Vstup: Zadání intervalu

4. Dosazení intervalu do f(x), ověření intervalu - rozhodování

nevyhovuje: Skok na bod 3. (opětovné zadání intervalu)

5. Deklarace proměnné n = -1 (stavová proměnná - počet iterací)

6. Vstup: Zadání požadované chyby

7. Cyklus Loop

8. Zvýšení proměnné n

9. Výpočet chyby (uložen do c)

10. Výpočet středu intervalu (uložen do s)

11. Zobrazení mezivýsledku

12. Ověření chyby

vyhovuje požadované chybě: skok na návěstí Ende (bod 16.)

13. Výpočet funkčních hodnot pro krajní bod a a střed intervalu

14. Ověření součinu f(a)*f(s)

f(a)*f(b) < 0 - střed s ulož do proměnné b (s›b)

f(a)*f(s) > 0 - střed s ulož do proměnné a (s›a)

f(a)*f(s)=0 - otestuj f(s)=0. Pokud platí, střed s ulož do a (s›a), skok na návěstí Ende (bod 16.)

15. Konec cyklu Loop (skok na začátek - bod 7.)

16. Návěstí Ende

17. Výstup: Zobrazení výsledku

18. Vymazání funkce f(x)

19. Zpětné nastavení kalkulačky

20. Konec programu

Zdrojový text programu v TI-Basicu

Program není příliš optimalizován. Prostředí pro výpočty je na začátku nastaveno na Approximate, protože se jedná o numerické výpočty. Také je deklarována skupina lokálních proměnných, ze kterých je vynechána proměnná f, ve které je uložena zkoumaná funkce. Učinil jsem tak záměrně, abych mohl porovnat výsledky s výsledkem výpočtu pomocí vestavěných funkcí (například zeros). Použití cyklu Loop se jeví nejvhodnější v případech, kdy není jisté, kolikrát cyklus proběhne. Zobrazování mezivýsledků je čistě demonstrační záležitostí. Řízení programu je provedeno pomocí příkazu If. Program je také k dispozici ke stažení

bisekce()
Prgm
Local a,b,c,d,n,fa,fb,fs,k,z
setMode("Exponential Format","NORMAL")
setMode("Exact/Approx","APPROXIMATE")

Lbl start
ClrIO
Input "Dolni mez",a
Input "Horni mez",b
f|x=a-›fa
f|x=b-›fb
fa*fb-›z
If b < a or z > 0 Then

	Disp "Spatne zvoleny interval"
	Pause
	Goto start
EndIf
-1-›n

Input "Zadej chybu",d
Loop
	n+1-›n
	ClrIO
	abs(b-a)/2^n-›c
	(a+b)/2-›s
	Disp "a,b,s=",{a,b,s}
	Disp "chyba",c
	Pause
	ClrIO
	If c < d Then
		Goto ende
	EndIf
	f|x=a-›fa
	f|x=s-›fs
	If fa*fs > 0 Then
		s-›a
	EndIf
	If fa*fs < 0 Then
		s-›b
	EndIf
	If fa*fb=0 Then
		If fs=0 Then
			s-›a
	EndIf
	Goto ende
	EndIf
EndLoop

Lbl ende
Disp "Vysledek:"

Disp string(a)&" +/- "&string(d)
setMode("Exact/Approx","AUTO")
EndPrgm