2.1 Aritmetika s pohyblivou řádovou čárkou

Kapitola 2

Aritmetika počítače V této kapitole se budeme zabývat vlivem zaokrouhlovacích chyb, které vznikají při numerických výpočtech prováděných na počítači v aritmetice s konečnou přesností. Bude nás zajímat, zda je algoritmus stabilní vůči zaokrouhlovacím chybám, tj. zda je výsledek výpočtu „dostatečně přesná aproximace řešení. Nejprve popíšeme vznik zaokrouhlovacích chyb a jejich šíření při provádění elementárních aritmetických operací.

2.1

{chp:2}

Aritmetika s pohyblivou řádovou čárkou

Číslo je v počítači zobrazeno jako posloupnost bitů (každý s číselným obsahem 0 nebo 1) konečné délky. Tato délka je pevně stanovena (např. 16, 32, 64 či 128 bitů). Počítač většinou umožňuje několik typů zobrazení čísel, jimž odpovídá několik velikostí paměťových míst. Nás především zajímá, jak jsou v počítači zobrazena reálná čísla. Je zřejmé, že při zvoleném typu zobrazení a předepsané délce paměťového místa je možno v počítači zobrazit pouze konečný počet čísel. Proto často hovoříme o konečné aritmetice či aritmetice s končenou přesností. Množina reálných čísel R je v počítači reprezentována svojí konečnou podmnožinou F , F ⊂ Q ⊂ R, kterou nazýváme soustavou čísel s pohyblivou řádovou čárkou (floating point number system), kde Q je množinou racionálních čísel. Prvky F lze zapsat ve tvaru (2.1) {cvut_2.1} y = ±m × β e−t , kde celé číslo β (obvykle β = 2) je nazýváno základem, celé číslo t určuje přesnost, celé číslo m pohybující se v rozsahu 0 ≤ m < β t − 1 je nazýváno mantisou a celočíselný parametr e exponentem. Množina F je plně určena parametry β, t a horní, resp. dolní mezí celočíselného exponentu emin ≤ e ≤ emax . Vztah 2.1 můžeme přepsat do názornějšího tvaru  d1 d2 dt e + 2 + . . . + t = ±β e × (0.d1 d2 . . . dt )β , y = ±β × (2.2) {cvut_2.2} β β β kde každá číslice di leží v intervalu 0 ≤ di ≤ β − 1 a (0.d1 d2 . . . dt )β představuje číslo zapsané v číselné soustavě se základem β. Je výhodné uvažovat m ≥ β t−1 pro y = 0, pak zřejmě d1 = 0 (nenulová první číslice mantisy). Systém dodržující tuto konvenci nazýváme normalizovaný. I když se v minulosti používaly různé základy β (do dnešní doby jsou rozšířeny základy β = 2 a β = 16) a stále se můžeme setkat s rozličnými hodnotami t, emin a emax , vývoj spěje k všeobecnému uznání tzv. IEEE standardní aritmetiky. Stučně ji vyložíme a další vlastnosti aritmetiky s pohyblivou řádovou čárkou budeme popisovat na tomto příkladu. 15

KAPITOLA 2. ARITMETIKA POČÍTAČE

16

IEEE aritmetika používá β = 2 a rozlišuje dva základní formáty čísel v pohyblivé řádové čárce: čísla s jednoduchou (single) a dvojitou přesností (double). V prvém případě je k uložení čísla použito 32, ve druhém 64 bitů. Uložení jednotlivých parametrů je patrné ze schematu na obrázku 2.1. V případě jednoduché přsnosti je na exponent vyhrazeno 8 bitů, do kterých je možno uložit celé číslo v rozmezí 0 až 255. Řetězce (0000 0000)2 = 0 či (1111 1111)2 = 255 mají však speciální význam, který bude popsán níže. Zbylá čísla 1 až 254 určují hodnotu veličiny e + 126, tj. hodnota exponentu e se pohybuje v rozmezí emin = −125 ≤ e ≤ 128 = emax . Na mantisu je vyhrazeno zbývajících 23 bitů, přičemž se standardně využívá normalizace, cifra d1 = 1 se přitom nezapisuje. S jednoduchou přesností uložené nenulové číslo v pohyblivé řádové čárce můžeme tedy zapsat ve tvaru y = ±2(g1 ...g8 )2 −126 × (0.1d2 d3 . . . d24 )2 ,

{cvut_2.3}

z čehož vyplývá

(2.3)

2−126 ≤ |y| ≤ (1 − 2−24 ) × 2128 ∼ 1038 .

Jednoduchá přesnost ±

0

exponent (g1 ...g8 )

mantisa (d2 ...d24 )

8

31

Dvojitá přesnost ±

0 {pic:sin_dou}

exponent (g1 ...g11 )

mantisa (d2 ...d53 )

11

63

Obrázek 2.1: Uložení parametrů IEEE standardni aritmetiky pro jednoduchou a dvojitou přesnost.

Čísla v pohyblivé řádové čárce nejsou vzhledem k R rovnoměrně rozložena (viz cvičení 2.1.1). Mají-li dvě čísla y1 , y2 ve vyjádření (2.3) shodný exponent e a jedná-li se o dvě po sobě jdoucí čísla množiny F , y2 > y1 , pak y2 − y1 = 2e × 2−24 . Rozložení čísel množiny F je charakterizováno pomocí strojové přesnosti M , což je vzdálenost čísla 1.0 od nejbližšího vyššího čísla v pohyblivé řádové čárce. Zřejmě platí M = 21 × 2−24 = 2−23 . Snadno lze ukázat, že vzdálenost libovolného normalizovaného čísla x ∈ F od svých sousedů je nejméně M |x|/2 a nejvýše M |x| (cvičení 2.1.2). Pokud by IEEE aritmetika pracovala pouze s množinou F , která obsahuje jen normalizovaná čísla popsaná (2.3), došlo by k nepříjemnému jevu. Zatímco čísla blízká 2−126 zprava jsou aproximována s chybou odpovídající počtu bitů mantisy, nejbližší číslo menší než 2−126 definované (2.3) je −2−126 , vzniká tak výrazná mezera v okolí nuly a zřejmě 0 ∈ F. K odstranění tohoto nedostatku množiny F

2.1. ARITMETIKA S POHYBLIVOU ŘÁDOVOU ČÁRKOU

17

přidává norma IEEE navíc množinu tzv. subnormálních čísel. To jsou nenulová nenormalizovaná čísla s exponentem (0000 0000)2 = 0, definovaná vztahem y = ±m × β emin −t , 0 < m < β t−1 , neboli

y = ±m × 2−149 , 0 < m < 223 .

Je-li řetězec mantisy i exponentu nulový (tj. znaménkový bit je 0 nebo 1, všechny ostatní bity jsou nulové) dostaneme reprezentaci čísel ±0 (+0 a −0 mají odlišnou reprezentaci, samozřejmě je zajištěno, že platí +0 = −0). Je-li řetězec exponentu roven (1111 1111)2 = 255 a mantisa je nulová, pak zobrazené číslo je definováno jako ±Inf (±∞). Je-li řetězec exponentu roven (1111 1111)2 = 255 a řetězec mantisy je nenulový (jeho hodnota je libovolná nenulová), pak je obsah interpretován jako NaN (Not a Number). Shrnutí je uvedeno v tabulce 2.1. řetězec exponentu (0000 0000)2 = (0)10 (0000 0001)2 = (1)10 .. .

numerická hodnota uloženého čísla ±(0.0d2 d3 . . . d24 )2 × 2−125 ±(0.1d2 d3 . . . d24 )2 × 2−125 .. .

(0111 1110)2 = (126)10 (0111 1111)2 = (127)10 .. .

±(0.1d2 d3 . . . d24 )2 × 20 ±(0.1d2 d3 . . . d24 )2 × 21 .. .

(1111 1110)2 = (254)10 (1111 1111)2 = (255)10

±(0.1d2 d3 . . . d24 )2 × 2128 ±∞ pokud d2 = d3 = . . . = d24 = 0; jinak NaN

Tabulka 2.1: IEEE jednoduchá přesnost.

{tab:single}

Zbývají dvě otázky. Jaká je přesnost zobrazení reálného čísla v množině F a jaká je přesnost provádění elementárních aritmetických operací +, −, ∗, /. Přesnost aproximace je charakterizována zaokrouhlovací jednotkou u = (1/2)β 1−t = (1/2)2−23 = 2−24 . Důkaz následující věty (která je formulována obecně) ponecháme na cvičení. Věta 2.1.1. Nechť x ∈ R leží mezi nejmenším a největším číslem množiny F . Označíme-li f l : R → F zobrazení z R do F , pak platí f l(x) = x(1 + δ), |δ| < u,

(2.4)

kde u je zaokrouhlovací jednotka.

{cvut_2.4} {veta_cvut_2.1}

Aritmetické operace +, −, ∗, / se v IEEE aritmetice v obvyklých případech provádějí tak, jako kdyby byly nejprve provedeny přesně (s nekonečnou přesností) a pak je výsledek zaokrouhlen na nejbližší číslo z F (v případě nerozhodnosti se zaokrouhluje dolů). Jsou-li x, y ∈ F, pak platí f l(x ± y) = (x ± y)(1 + δ1 ), |δ1 | ≤ u, f l(x ∗ y) = (x ∗ y)(1 + δ2 ), |δ2 | ≤ u,

(2.5) {cvut_2.5}

f l(x/y) = (x/y)(1 + δ3 ), |δ3 | ≤ u. Analogický vztah se obvykle přepokládá i pro operaci odmocnění. Nastane-li výjimečný případ, je výsledek generován podle tabulky 2.2. Přetečením rozumíme případ kdy je přesný výsledek operace v absolutní hodnotě větší než největší číslo z F . Podtečením rozumíme případ, kdy je přesný výsledek operace v absolutní hodnotě menší, než nejmenší kladné normalizované číslo.

KAPITOLA 2. ARITMETIKA POČÍTAČE

18 typ výjimky nedefinované operace přetečení dělení nenulového čísla nulou podtečení

příklad √ 0/0, 0 × ∞, −1

výsledek NaN ±∞ ±∞ subnormální čísla

Tabulka 2.2: Výjimky v IEEE aritmetice.

{tab:exceptns}

Vlastnosti aritmetiky s pohyblivou řádovou čárkou jsme vyložili na příkladu IEEE aritmetiky s jednoduchou přesností. Je zřejmé, jak postupovat při odvození charakteristik aritmetiky založené na jiné hodnotě parametrů. Pro doplnění uvádíme v tabulce 2.3 porovnání IEEE aritmetiky s jednoduchou a dvojitou přesností. přesnost jednoduchá dvojitá

počet bitů celkem mantisa 32 23(+1) 64 52(+1)

exp. 8 11

zaokrouhlovací jednotka u 2−24 ∼ 5.96 × 10−8 2−53 ∼ 1.11 × 10−16

rozsah 10±38 10±308

Tabulka 2.3: IEEE aritmetika s pohyblivou řádovou čárkou, jednoduchá a dvojitá přesnost. {tab:sin_dou}

Cvičení – konečná aritmetika {trivaritmetic} {trivaritmetic2}

{trivaritmetic3}

Doporučená literatura: [11], [16], [4]. Příklad 2.1.1. Vypočtěte a graficky znázorněte na číselné ose prvky množiny čísel s pohyblivou řádovou čárkou pro β = 2, t = 3, emin = −1 a emax = 3. Příklad 2.1.2. Nechť je množina F ⊂ R konečnou množinou reprezentující množinu reálných čísel v počítači, množinou normalizovaných čísel. Ukažte, že vzdálenost Δx libovolného normalizovaného čísla x z množiny F od jeho nejbližšího souseda je M |x| ≤ Δx ≤ M |x|. β −1 M |x| = 2 Příklad 2.1.3. Dokažte větu 2.1.1. Nechť x ∈ R leží mezi nejmenším a největším číslem množiny F . Označíme-li f l : R → F zobrazení z R do F , pak platí f l(x) ≤ x(1 + δ), |δ| < u

{trivaritmetic4}

{vyrokyvIEEE}

kde u = M /2 se nazývá zaokrouhlovací jednotka. Příklad 2.1.4. Odvoďte parametry IEEE aritmetiky s pohyblivou řádovou čárkou v dvojité přesnosti (64 je celkový počet bitů, z toho 52(+1) bitů tvoří matnisu, 11 bitů tvoří exponent, viz tabulku 2.3). Příklad 2.1.5. Který z následujících výroků je pravdivý v IEEE aritmetice, předpokládáme-li, že a, b a c jsou normalizovaná čísla v pohyblivé řádové čárce a že nenastane žádná výjimečná situace (například přetečení) (a) f l(a op b) = f l(b op a), kde op = +, ×

2.2. ZAOKROUHLOVACÍ CHYBY V ARITMETICE S KONEČNOU PŘESNOSTÍ19 (b) f l(b − a) = −f l(a − b) (c) f l(a + a) = f l(2 · a) (d) f l(0.5 · a) = f l(a/2) (e) f l((a + b) + c) = f l(a + (b + c)) (f) je-li a ≤ b, pak a ≤ f l((a + b)/2) ≤ b Případnou nepravdivost výroků zkuste ověřit nalezením příkladu. Výhodně k tomu můžete využít jednoduché aritmetiky z příkladu 2.1.1.

2.2

Zaokrouhlovací chyby v aritmetice s konečnou přesností

Při povrchním pohledu na vztah (2.4) a (2.5) by se mohlo zdát, že zaokrouhlovací chyby jsou velmi malé a jejich vliv při provádění numerických výpočtů nebude velký (snad s výjimkou velkého počtu operací s nějakými extrémními čísly). Ukážeme na několika příkladech, že tento ukvapený závěr je zcela mylný. Prvním příkladem je tzv. krácení platných cifer (cancellation), které nastává, odečítáme-li dvě téměř shodná čísla. Příklad. Uvažujme funkci f (x) = (1 − cos(x))/x2 (převzato z [11]). Pro x = 1.2 × 10−5 je hodnota cos(x) zaokrouhlená na 10 desetinných míst rovna c = 0.9999999999, vyčíslením hodnoty f (1.2 × 10−5 ) dostaneme (1 − c)/x2 = 10−10 /(1.44 × 10−10 ) = 0.6944 . . . , což je úplně špatně, neboť 0 ≤ f (x) ≤ 0.5 pro x = 0. Vidíme, že přestože hodnota cos(x) byla aproximována s přesností na 10 desetinných míst, výsledek výpočtu f (x) neaproximuje správnou hodnotu ani s přesností jednoho desetinného místa! Je důležité si uvědomit, že problém není způsoben vlastním odečtením (1 − c), to bylo provedeno přesně. Problém spočívá v tom, že sama hodnota c byla určena nepřesně a výsledek přesného výpočtu (1 − c) je díky krácení platných cifer stejného řádu, jako je chyba obsažená v c. Tím se významnost nepatrné chyby posunula o 10 řádů a katastrofálně ovlivnila celý další výpočet, byť byl proveden sebepřesněji (často se proto hovoří v této souvislosti o tzv. katastrofickém krácení). Pokusíme se krácení popsat pomocí vztahů (2.4) a (2.5). Nechť x ˆ a yˆ jsou dvě čísla zažížená jistou chybou, tj x ˆ = x(1 + δx), yˆ = y(1 + δy). Předpokládejme, že chyby δx, resp. δy jsou malé vzhledem k velikosti x, resp. y, může jít o chyby způsobené předcházejícím výpočtem nebo třeba o zaokrouhlovací chyby při uložení dat do počítače (pak x ˆ = f l(x), yˆ = f l(x) a |δx| ≤ u, |δy| ≤ u). Proveďme přesný součet čísel x ˆ a yˆ (čísla mohou mít opačná znaménka, příklad zahrnuje i odečítání) sˆ = x ˆ + yˆ = x(1 + δx) + y(1 + δy) = = x + y + x δx + y δy = = (x + y)(1 + δs),

KAPITOLA 2. ARITMETIKA POČÍTAČE

20 kde δs =

y x δx + δy. x+y x+y

Je zřejmé, že i když hodnoty δx a δy jsou malé, není zaručeno, že hodnota δs bude rovněž malá. Pokud bude x  (x + y) a zároveň δx = 0, respektive y  (x + y) a zároveň δy = 0, bude chyba δs relativně velká, δs  δx, respektive δs  δy. Znovu vidíme, že krácení není nebezpečné samo o sobě (dojde-li ke krácení při odečtení dvou přesných hodnot, k žádné ztrátě přesnosti nedojde), ale je nebezpečné tím, že zesiluje vliv předchozích chyb obsažených v datech. Druhý příklad ukazuje, že i bez krácení popsaného výše, může dojít při provedení jednoduchého výpočtu k velké chybě. Příklad. Předpokládejme, že chceme nalézt dobrou numerickou aproximaci hodnoty e s použitím vztahu e = limn→∞ (1 + 1/n)n , kde limitu nahradíme prostým výpočtem hodnoty f (n) = (1 + 1/n)n pro dostatečně velké n. Pro n = 10 dostaneme v případě výpočtu v IEEE aritmetice s jednoduchou přesností lepší aproximaci čísla e než pro n = 107 ! Příčina je následující. Sčítáme-li (1 + 1/n) pro n  1, obsahuje výsledek součtu stále méně a méně informace o čísle n (neboť 1  1/n). I když provedeme následné umocnění přesně, výsledek je zatížen velkou chybou, viz cvičení, příklad 2.2.2. Posledním příkladem je sčítání řad s kladnými členy. Příklad. Z teorie Fourierových řad je známo, že ∞

n−2 = π 2 /6.

n=1

Předpokládejme, že tuto identitu neznáme a chceme vypočítat součet řady numericky sčítáním (. . . (((1 + 2−2 ) + 3−2 ) + 4−2 ) + . . . + m−2 ), kde m určíme jako nejmenší celé číslo, jehož zahrnutí do výpočtu již nezmění vypočtený součet. Výsledek výpočtu bude překvapivě nepřesný (viz cvičení, příklad 2.2.3). Příčina je opět zřejmá: řada konverguje velmi pomalu a náš výpočet je prováděn tak, že hodnota přičítaných prvků se stále zmenšuje a hodnota částečných k −2 n takový, součtů se stále zvětšuje. Pro jisté k je vypočtený částečný součet n=1 ∞ −2 −2 že přičtení (k + 1) nezmění jeho hodnotu. Zbytek n=k+1 n je však stále příliš velký. Jak překonat popstanou obtíž? První nápad může být změnit pořadí sčítání (sčítat od nejmenšího prvku k největšímu). Problém ovšem je, že nevíme, kterým prvkem začít. Navíc, uspořádání sčítanců je obecně drahá operace a nelze ji v praktických výpočtech použít. Univerzálním řešením je použití speciálních technik zvyšujících přesnost (samozřejmě na úkor rychlosti). Zvídavého čtenáře odkazujeme na [11, Chapter 4]. Jiným řešením může být použití vhodné identity a řady konvergující podstatně rychleji, viz cvičení, příklad 2.2.6. V každém případě je vhodné zamyslet se nad konvergencí sčítané řady. Odstrašujícím případem budiž všem eventuální pokus nalézt ∞ výše popsaným postupem součet řady n=1 n−1 .

2.2. ZAOKROUHLOVACÍ CHYBY V ARITMETICE S KONEČNOU PŘESNOSTÍ21

Cvičení – rušení platných číslic, součty řad Doporučená literatura: [3], [4]. Příklad 2.2.1. Ukažte, jak je potřeba přepsat uvedené výrazy, aby byl omezen vliv krácení platných cifer: √ (a) pro x ∼ 0 výraz x + 1 − 1, (b) pro x ∼ y výraz sin(x) − sin(y), (c) pro x ∼ y výraz x2 − y 2 , (d) pro x ∼ 0 výraz (1 − cos(x))/sin(x). Příklad 2.2.2. Vypočtěte hodnotu výrazu (1 + 1/n)n pro n = 101 , 102 , . . . , 107 a srovnejte výsledky s číslem e. ∞ Příklad 2.2.3. Vypočtěte aproximaci součtu nekonečné řady n=1 1/n2 porovnejte s hodnotou π 2 /6. Určete chybu a udejte, kolik členů řady jste sčítali. ∞ 2 Příklad 2.2.4. Vypočtěte aproximaci součtu nekonečné řady n=1 1/(n + 1). Součet řady aproximujte k-tým částečným součtem. Výpočet zastavte při splnění podmínky  k  k+1 

1

1 fl = fl . n2 + 1 n2 + 1 n=1 n=1 Určete počet členů k řady, které jste sčítali. ∞ 2 Příklad 2.2.5. Vypočtěte aproximaci součtu nekonečné řady n=1 1/(n + 1) −6 s přesností větší než 10 . Při určení počtu členů řady k použijte nerovnost ∞

n=k+1

1 ≤ n2 + 1



∞

k

{priklad_pi26} {sum_naiv}

{sum_presnost}

π dx = arctan(∞) − arctan(k) = − arctan(k) < 10−6 . x2 + 1 2

Sčítejte od nejmenších členů k největším (tedy od k do 1). Proč je to výhodnější? ∞ 2 Příklad 2.2.6. Vypočtěte aproximaci součtu nekonečné řady n=1 1/(n + 1) −6 s přesností větší než 10 . Nesčítejte však původní řadu, ale využijte následujících identit ∞ ∞ ∞

1 π2 1 π4 1 π6 , , . = = = 2 4 6 n 6 n=1 n 90 n=1 n 945 n=1 Opět sčítejte od nejmenších členů k největším. Návod k řešení. Zřejmě platí následující vztahy: x≡

{priklad_e}

∞

1 , 2+1 n n=1

x1 ≡

∞ ∞

1 1 π2 1 = = , − 2 2 2 2 6 − x n=1 n n + 1 n=1 n (n + 1)

x2 ≡

∞ ∞

π2 1 1 1 = − = , 4 4 + n2 4 (n2 + 1) 90 − x1 n n n n=1 n=1

∞ ∞

1 1 π2 1 . = − 6 = x3 ≡ 6 4 6 (n2 + 1) 945 − x2 n n + n n n=1 n=1

{sum_identity}

KAPITOLA 2. ARITMETIKA POČÍTAČE

22

Výpočet proveďte nejpve s využitím pouze první identity, poté s využitím prvních dvou identit a nakonec s využitím všech tří identit. Výsledek určete nejprve „naivním způsobem (jako v příkladu 2.2.4), dále sčítejte každou řadu od nejmenších členů k největším se zaručenou přesností (viz příklad 2.2.5). Při určení počtu sčítanců k v jednotlivých případech použijte odhady pomocí vhodných integrálů. Platí následující vztahy  dx = arctan(x), x2 + 1   dx 1 dx = − , n ≥ 2. − n 2 n−1 n−2 x (x + 1) (n − 1)x x (x2 + 1) Dále je vhodné využít nerovnosti  b a

dx ≤ xn (x2 + 1)

 a

b

dx . xn+2

Porovnejte výsledky z příkladů 2.2.4, 2.2.5, 2.2.6 (jednotlivé součty a počty sčítanců, ev. dobu, kterou počítač k vyčíslení součtů potřeboval počítal). ∞ 2 Příklad 2.2.7. Zkuste vypočítat aproximaci součtu řady n=1 1/(n + 1) , nebo nějaké jiné vámi zvolené nekonečné řady. Aproximaci součtu vypočtěte opět několika různými způsoby, tak jako u předchozích příkladů. Příklad ∞2.2.8. Zkuste „naivním způsobem vypočítat aproximaci součtu nekonečné řady n=1 1/n. Jaká je asi přesnost takto získaného výsledku? ∞ Návod k řešení. Pro řadu n=1 1/n je podmínka f l(sk ) = f l(sk+1 ) splněna až pro velmi vysoké k. Zamyslete se nad konvergencí této řady a zkuste raději „naivním způsobem ∞ aproximovat součet jiné nekonečné řady s obdobnými vlastnostmi, například n=1 1/(n − 0.9999999995).

Cvičení – příklady v MATLABu Doporučená literatura: [19], [4]. Příklad 2.2.9. Máme následující program v MATLABu >> a = 1; >> b = 1; >> while a+b ~= a; b=b/2; end; Promněnná b je v cyklu půlená, dokud je součet a + b různý od a. Budeme-li výpočet provádět na množině reálných čísel, program skončí v nekonečném cyklu. Náš program však skončí po konečném počtu kroků. Proveďtě výpočet. Určete hodnotu b, při které program skončí. O jakou hodnotu se jedná? Návod k řešení. Jedná se o zaokrouhlovací jednotku u = M /2 = 1.1102 × 10−16 . Relativní strojovou přesnost M v MATLABu zjistíme příkazem eps.

2.2. ZAOKROUHLOVACÍ CHYBY V ARITMETICE S KONEČNOU PŘESNOSTÍ23 Příklad 2.2.10. Nyní budeme krátkým výpočtem dokumentovat neplatnost zákona asociativity sčítání v aritmetice čísel s pohyblivou řádovou čárkou. Máme následující tři čísla >> a = 1.0e+308; >> b = 1.1e+308; >> c = -1.001e+308; Součet čísel a + b + c spočítáme v MATLABu dvěma různýmy způsoby. Dostaneme následující výsledky. >> a+(b+c) ans = 1.0990e+308

>> (a+b)+c ans = Inf

Vysvětlete. Jaký je rozdíl mezi tímto příkladem a příkladem 2.1.5 (e)? Návod k řešení. V druhém výpočtu došlo k výjimce přetečení. Příklad 2.2.11. Mějme výraz (1 + x) − 1 , x hodnota tohoto výrazu je zřejmě rovna 1 pro libovolné x = 0. Pro pevně zvolené x jsme výraz vyčíslili v MATLABu. Následuje výpis programu >> x = 1.0e-15; >> ((1+x)-1)/x ans = 1.1102 Kde nastala chyba? Návod k řešení. Rušení platných číslic – cancelation.

Cvičení – polynomy Doporučená literatura: [11], [4], [3]. Příklad 2.2.12. Metodou bisekce nalezněte kořen(y) polynomu p(x) = x9 − 18x8 + 144x7 − 672x6 + 2016x5 − 4032x4 + 5376x3 − 4608x2 + 2304x− 512 ležící v intervalu [1.8, 2.2]. Jeden z kořenů polynomu by měl ležet blízko hodnoty 2. Zkuste ověřit vámi získaný výsledek aplikováním metody jiné intevaly blízké intervalu původnímu, např. [1.9, 2.2], [1.7, 2.2]. Vykreslete graf tohoto polynomu v okolí nalezených kořenů. K vyčíslení hodnoty polynomu d

ai xi p(x) = i=0

použijte Hornerova pravidla. V MATLABu určíme hodnotu polynomu p(x) v bodu x programem

{polynom}

KAPITOLA 2. ARITMETIKA POČÍTAČE

24 >> >> >> >>

coeff = [1 -18 144 -672 2016 -4032 5376 -4608 2304 -512]; n = 10; p = 0; for i = 1:n, p = x*p + coeff(i); end

MATLAB obsahuje funkci polyval(), která implementuje Hornerovo pravidlo, odpovídající program pak vypadá >> coeff = [1 -18 144 -672 2016 -4032 5376 -4608 2304 -512]; >> p = polyval(a,x);

{polynom_roots}

Konfrontujte graf s výsledky získanými metodou bisekce, vysvětlete. Příklad 2.2.13. Pro výpočet kořenů polynomu je v MATLABu přímo implementovaná funkce roots(). Zkuste kořeny polynomu z předchozího příkladu nalézt pomocí této funkce >> coeff = [1 -18 144 -672 2016 -4032 5376 -4608 2304 -512]; >> r = roots(coeff); Porovnejte výsledky s výsledky, které jste získali metodou bisekce v předchozím příkladu. Pokuste se vysvětlit. Příklad 2.2.14. Zopakujte postupy z předchozích dvou příkladů pro jednodušší polynom q(x) = x7 − 7x6 + 21x5 − 35x4 + 35x3 − 21x2 + 7x − 1.

{priklad_raclomfce}

Při hledání kořenů pomocí bisekce předpokládejte, že jeden z kořenů polynomu leží blízko hodnoty 1. Sami navrhněte polynom r(x), u něhož znáte alespoň přibližně hodnotu minimálně jednoho z jeho kořenů, a aplikujte výše zmíněné postupy. Srovnejte chování všech polynomů p(x), q(x) a r(x). Čím jsou polynomy p(x) a q(x), popřípadě i r(x), charakteristické? Příklad 2.2.15. Vykreslete graf racionální lomené funkce 4x4 − 59x3 + 324x2 − 751x + 622 . x4 − 14x3 + 72x2 − 151x + 112 Hodnoty polynomů čitatele a jmenovatele počítejte pomocí Hornerova pravidla. Graf funkce vykreslete pro hodnoty ρ(x) =

x = a + (k − 1) × 2−26 , kde a = 1.606 a k = 1, . . . , 361 (příklad převzat z [11]). dρ Spočtěte derivaci π(x) = dx funkce ρ(x) a její numerickou aproximaci (z hodnot funkce ρ(x)). Vykreslete grafy a porovnejte. Výpočty provádějte vždy v jednoduché a dvojité přesnosti. Zhodnoťte výsledky. Návod k řešení. Pozor, některé verze MATLABu neumí operovat s čísly v jednoduché přesností. Při vyčíslování funkce v jednoduché přesnosti je tedy potřeba provádět konverzi příslušných mezivýsledků pomocí funkcí single() a double().