´ Uvod do numerick´ e matematiky 1
Pˇ redmˇ et numerick´ e matematiky
Numerick´ a matematika je vˇeda, kter´a se zab´yv´a ˇreˇsen´ım matematicky formulovan´ych ´uloh pomoc´ı logick´ych operac´ı a aritmetick´ych operac´ı s ˇc´ısly o koneˇcn´e d´elce. 2 typy matematicky formulovan´ ych u ´ loh • numericky formulovan´e ´ulohy - jednoznaˇcn´y funkˇcn´ı vztah mezi koneˇcn´ym poˇctem vstupn´ıch a v´ystupn´ıch dat, jedn´a se obvykle o algebraick´e u ´lohy, nˇekdy je moˇzno nal´ezt teoretick´e ˇreˇsen´ı ´ulohy pomoc´ı koneˇcn´e posloupnosti aritmetick´ych a logick´ych operac´ı, jindy ne (lze nal´ezt pouze pˇribliˇzn´e ˇreˇsen´ı) • u ´lohy, kter´e nejsou numericky formulovan´e - obvykle ´ulohy matematick´e anal´yzy, ve kter´ych je obsaˇzen nekoneˇcnˇe kr´atk´y krok (napˇr. v´ypoˇcet derivace, integr´alu, ˇreˇsen´ı diferenci´aln´ı rovnice); tyto ´ulohy je tˇreba nejdˇr´ıve nˇejak´ym zp˚ usobem pˇrev´est na numerick´e ´ulohy ´lohy nebo jej´ı Numerickou metodou rozum´ıme postup v´ypoˇctu numerick´e u pˇrevod na ´ulohu jednoduˇsˇs´ı ˇci postup, kter´y nahrazuje matematickou u ´lohu ´ulohou numerickou. Algoritmem rozum´ıme realizaci numerick´e metody, tj. konkr´etn´ı koneˇcnou posloupnost operac´ı, kter´a s poˇzadovanou pˇresnost´ı pˇrevede vstupn´ı data na v´ysledn´e hodnoty. Algoritmus lze programovat na poˇc´ıtaˇci.
1
2
Chyby - nepˇ resnosti pˇ ri ˇ reˇsen´ı u ´ loh
Zdroje chyb 1. Chyby vstupn´ıch dat (napˇr. chyby mˇeˇren´ı, chyby modelu reality) 2. Chyby metody (Truncation errors) - v d˚ usledku pˇreveden´ı matematick´e ´ulohy na numerickou 3. Zaokrouhlovac´ı chyby (Roundoff errors) - v d˚ usledku zaokrouhlov´an´ı pˇri v´ypoˇctech s ˇc´ısly o koneˇcn´e d´elce Definice chyb x - pˇresn´a hodnota
x˜ - pˇribliˇzn´a hodnota
A(x) = |˜ x − x| ≤ a(x) A(x) - absolutn´ı chyba R(x) =
a(x) - odhad absolutn´ı chyby
A(x) ≤ r(x) |x|
R(x) - relativn´ı chyba
−→ ∼
r(x) ≃
a(x) .... pokud r(x) ≪ 1 |˜ x|
r(x) - odhad relativn´ı chyby
Intervalov´y odhad x x˜ − a(x) ≤ x ≤ x˜ + a(x) −→ x = x˜ ± a(x) x = x˜ · (1 ± r(x)) Relativn´ı chyba ⇐⇒ Poˇ cet platn´ ych ˇ c´ıslic Necht’ x 6= 0, ˇc´ıslo x zap´ıˇseme pomoc´ı ˇc´ıslic x1 , x2, . . . , xn, necht’ m ≤ n x = ±0.x1x2 . . . xm−1xm xm+1 . . . xn · 10p
necht’ x1, . . . , xm jsou platn´e cifry, x1 6= 0. Pak r(x) < 5.10−m .....( pˇresnˇeji
r(x) < 5.10−m/x1)
Pˇresnost v´ypoˇctu je obvykle d´ana relativn´ı chybou. 2
3
Chyby metody (aproximace) • Pˇri v´ypoˇctech derivace, integr´alu apod. nahrazujeme nekoneˇcnˇe kr´atk´y krok dx koneˇcn´ym krokem h • Tento typ chyby nijak nesouvis´ı se zaokrouhlov´an´ım • 1 veliˇcinu lze aproximovat mnoha r˚ uzn´ymi zp˚ usoby
ˇ ad metody - Je-li chyba δy veliˇciny y ´umˇern´a δy ∼ hα ∼ O(hα ) • R´ pak α naz´yv´ame ˇr´adem metody
Uk´ azky r˚ uzn´ ych zp˚ usob˚ u aproximace - odvozen´ı chyby metody (obvykle se pouˇ z´ıv´ a Taylor˚ uv rozvoj) derivace df y= dx x1
−→ y ≃
h f (x1 + h) − f (x1) = f ′(x1) + f ′′(x1) + · · · h 2
f ′′(x1) δy = h+ ··· 2
.....
metoda 1. ˇr´adu
df h2 ′′′ f (x1 + h/2) − f (x1 − h/2) ′ y= = f (x1) + f (x1) + · · · −→ y ≃ dx x1 h 24
f ′′′(x1) 2 h + ··· δy = 24
.....
metoda 2. ˇr´adu
3
ˇs´ıˇrka intervalu h = b−a N
integr´al - obd´eln´ıkov´a metoda y=
Z b
f (x)dx −→ y ≃ h
a
N X
f (xi) ..........
xi = a +
i=1
2i − 1 h 2
chyba v 1 podintervalu Z x +h/2 i xi −h/2
t2 ′′ ′ f (x) dx = f (xi) + t f (xi) + f (xi) + · · · dt = −h/2 2 h3 ′′ = hf (xi) + 0 + f (xi) 12 Z h/2
celkov´a chyba |δy| ≃
3 N h X ′′ f (xi) 12 i=1
|b − a| h3 N max |f ′′ (xi)| ≤ max |f ′′ (x)|h2 ∼ O(h2 ) ≤ i∈h1,N i x∈(a,b) 12 12
Chyba metody vyˇsˇs´ıho ˇr´adu kles´a rychleji pˇri zmenˇsov´an´ı kroku h. Pokud jsou metody jinak rovnocenn´e, vybereme metodu vyˇsˇs´ıho ˇr´adu. Znalost ˇr´adu metody umoˇzˇnuje • Odhad chyby • Zpˇresnˇen´ı v´ysledku Nejjednoduˇsˇs´ı zp˚ usob - spoˇctu odhady v´ysledku yh pro krok h a yh/2 pro krok h/2 Pˇr´ıklad pro metodu 1.ˇr´adu yh = y + ah + bh2 ! h h 2 yh/2 = y + a + b 2 2 Koeficienty a,b ˇcasto nemohu spoˇc´ıtat, ale pˇresto plat´ı δyh/2 ≃ yh − yh/2
......
(≃ ah/2 )
b 2 h = y + O(h2 ) 2 je chyba odhadnuta a ˇr´ad metody o 1 zv´yˇsen.
y˜ = 2yh/2 − yh = y − a tedy kombinac´ı yh a yh/2
4
4
Zaokrouhlovac´ı chyby
4.1
Reprezentace ˇ c´ısla v poˇ c´ıtaˇ ci
1. Cel´ aˇ c´ısla - pˇresn´e v´ypoˇcty, velmi omezen´y rozsah • INTEGER - 2 byty (INTEGER*2, v C++ short int) 16 bit˚ u = 216 ˇc´ısel h−32768, 32767i • LONGINT - 4 byty D(INTEGER*4,E v C++ int) 32 bity = 232 ˇc´ısel −231, 231 − 1 2. Re´ aln´ aˇ c´ısla - ˇc´ısla v pohybliv´e desetinn´e teˇcce - FLOATING POINT ≃ vˇedeck´y tvar ˇc´ısla ✏ ✓ ± 23 ✗✔ ✓
5.423687
±
✒ ✖✕ ✡ ✡ ✡ ✢
Znam´enko
✄ ✄✎
✄
Mantisa
✏
✒
.10
✑
✄ ✎✄
✄
✄ ✄
✑
Exponent
Mantisa a exponent jsou v poˇc´ıtaˇci bin´arn´ı ˇc´ısla, z´aklad u exponentu je 2. u na mantisu =⇒ pˇresnost ˇc´ısla D´ elka MANTISY - tj. poˇcet bit˚ pˇresnost ⇐⇒ poˇcet ˇc´ısel mezi 1 a 2 interval mezi ˇc´ısly mezi 1 a 2 je rovnomˇern´y - do pamˇeti se mohou ukl´adat ˇ ım v´ıce bit˚ jenom ˇc´ısla 1, 1 + ε, 1 + 2 ε, · · ·, 2 - ε. C´ u na mantisu, t´ım menˇs´ı ε =⇒ menˇs´ı chyby pˇri zaokrouhlov´an´ı (u meziv´ysledk˚ u je v registrech procesoru pˇresnost vyˇsˇs´ı). Byly zde uvedeny vˇsechny mantisy, pˇri zmˇen´ach zmˇen´ach exponent˚ u se exponent krok mezi ˇc´ısly zv´yˇs´ı ´umˇernˇe 2 (relativn´ı chyba ˇc´ısla se ale nemˇen´ı). u na exponent - urˇcuje rozsah D´ elka EXPONENTU - tj. poˇcet bit˚ Pozn. 1 exponent se zvl´aˇst’ mus´ı vyhradit pro 0, kter´a nem´a logaritmus, mantisy u tohoto exponentu lze vyuˇz´ıt pro vyznaˇcen´ı chyb (overflow, undefined), d´ale se m˚ uˇze vyuˇz´ıt pro ˇr´ıdkou s´ıt’ ˇc´ısel pod minimem k oˇsetˇren´ı podteˇcen´ı (underflow)
5
8 - 11 bit˚ u na exponent 7 8 bit˚ u na exponent napˇr. h−126, 127i − 22 ≃ 3.4 . 1038 8 9 bit˚ u na exponent - 22 ≃ 1.15 . 1077 9 10 bit˚ u na exponent - 22 ≃ 1.32 . 10154 10 11 bit˚ u na exponent - 22 ≃ 1.74 . 10308 Jednoduch´a pˇresnost = 4 byty C++ - float Fortran - Real = Real*4 norma IEEE urˇcuje zp˚ usob zaokrouhlov´an´ı, ne bity na mantisu a exponent ˇcasto 1 bit znam´enko + 8 bit˚ u exponent + 23 bit˚ u mantisa ⇒ ε ≃ 1, 2 . 10−7 u Dvojit´a pˇresnost = 8 byt˚ C++ - double Fortran - Double = Real*8 ˇcasto 1 bit znam´enko + 11 bit˚ u exponent + 52 bit mantisa ⇒ ε ≃ 2, 2 . 10−16 Dalˇs´ı typy Fortran - Complex, Complex*16, Real*16 TurboPascal - Real (6 byt˚ u, pˇresnost 1.5) TurboPascal - Extended (Real*10), Comp (Integer*8)
6
4.2
ˇ ıˇ S´ ren´ı chyb ve v´ ypoˇ ctech
Nebezpeˇcn´e jsou operace, kter´e mohou podstatnˇe zvˇetˇsit relativn´ı chybu !! • Sˇc´ıt´an´ı, odeˇc´ıt´an´ı a(x ± y) = a(x) + a(y)
−→
r(x ± y) =
a(x) + a(y) |x ± y|
ˇ Pokud v´ysledek mal´y =⇒ zvˇetˇs´ı se silnˇe relativn´ı chyba r !! Casto nemohu rozhodnout, zda v´ysledek je 0 nebo nen´ı !! Odeˇcteme-li napˇr. ˇc´ısla 1.32483726, 1.32483357 zn´am´a na 9 platn´ych ˇc´ıslic r(x) = r(y) ≃ 5 · 10−9 (pˇresnˇeji 3.8 · 10−9) dostaneme v´ysledek x − y = 0.00000369 s pˇresnost´ı na 3 platn´e ˇc´ıslice r(x−y) ≃ 5·10−3 (pˇresnˇeji r(x−y) = 1 · 10−3). Motivace v´yvoje ˇrady numerick´ych postup˚ u - snaha vyhnout se odeˇc´ıt´an´ı dvou pˇribliˇznˇe stejnˇe velk´ych ˇc´ısel. • N´asoben´ı, dˇelen´ı a(x · y) = |x|a(y) + |y|a(x) ! |x|a(y) + |y|a(x) x = a y y2
−→ −→
r (x · y) = r(x) + r(y) ! x = r(x) + r(y) r y
N´asoben´ı a dˇelen´ı nemohou podstatnˇe zvˇetˇsit zaokrouhlovac´ı chybu, nejsou tedy nebezpeˇcn´e. Pozn. Dˇelen´ı ˇc´ıslem 0 je hrub´a chyba - nejde o zaokrouhlovac´ı chybu.
7
4.3
Z´ avislost charakteru chyby na velikosti kroku h r
✻
✲ ✛
✲
dominuje chyba
dominuje chyba
zaokrouhlovac´ı
metody
h
Zaokrouhlovac´ı chyby pˇri mal´em kroku h vznikaj´ı z r˚ uzn´ych pˇr´ıˇcin - u derivace v d˚ usledku odeˇc´ıt´an´ı pˇribliˇznˇe stejn´ych ˇc´ısel, u integr´alu v d˚ usledku velk´eho poˇctu operac´ı Pozn. V obou pˇr´ıpadech pˇri stanoven´ı pˇr´ıliˇs kr´atk´eho kroku h m˚ uˇze doj´ıt i k hrub´ym chyb´am vzhledem k rozd´ılu rovn´emu 0 v d˚ usledku zaokrouhlen´ı nebo vzhledem k opakovan´emu prov´adˇen´ı souˇct˚ u, kter´e se pˇri dan´e numerick´e pˇresnosti neprojev´ı u + δu = u !
8
5
Korektnost a podm´ınˇ enost u ´ lohy
Korektnost u ´ lohy Definice: Necht’ ´ulohou je naj´ıt ˇreˇsen´ı ~y ∈ N (N je mnoˇzina moˇzn´ych ˇreˇsen´ı) pro zadan´y vektor ~x ∈ M (M je mnoˇzina vstupn´ıch dat). Pak u ´loha je korektn´ı pr´avˇe tehdy, jsou-li z´aroveˇn splnˇeny obˇe n´asleduj´ıc´ı podm´ınky 1. ∃ pr´avˇe jedno ˇreˇsen´ı ~y pro ∀ ~x ∈ M.
ˇ sen´ı spojitˇe z´avis´ı na vstupn´ıch datech, tj. jestliˇze pro ∀n z mnoˇziny 2. Reˇ pˇrirozen´ych ˇc´ısel je ~yn ˇreˇsen´ı pro vstupn´ı data ~xn , a jestliˇze ~y je ˇreˇsen´ı pro vstupn´ı data ~x, necht’ d´ale ρ je norma v mnoˇzinˇe vstupn´ıch dat a σ je norma v mnoˇzinˇe moˇzn´ych ˇreˇsen´ı, pak plat´ı ρ
σ
~xn → ~x ⇒ ~yn → ~y V praxi se ˇreˇs´ı i nekorektn´ı u ´lohy, ale 1. krok ˇreˇsen´ı spoˇc´ıv´a v nalezen´ı vhodn´eho zp˚ usobu, jak pˇrev´est u ´lohu na u ´lohu korektn´ı (napˇr. podm´ınkou na v´ysledek; interpretac´ı vstupn´ıch dat; vhodnou volbou normy v prostoru ˇreˇsen´ı apod.)
9
Podm´ınˇ enost u ´ lohy Definice: Podm´ınˇenost ´ulohy Cp je dan´a pomˇerem relativn´ı zmˇeny v´ysledku ku relativn´ı zmˇenˇe vstupn´ıch dat, tj. Cp =
kδyk kyk kδxk kxk
≈
r(y) r(x)
Pokud Cp mal´e (Cp ∼ 1), ˇr´ık´ame, ˇze ´uloha je dobˇre podm´ınˇen´a, pokud je Cp velk´e (Cp > 100), ´uloha je ˇspatnˇe podm´ınˇen´a. Pokud je pˇresnost pouˇzit´eho typu ˇc´ısel ε (r(x) = ε) , pak ´uloha s Cp > ε−1 nen´ı v r´amci dan´e pˇresnosti ˇreˇsiteln´a. ˇ Casto se pro ˇspatnˇe podm´ınˇen´e ´ulohy pouˇz´ıvaj´ı speci´aln´ı metody, kter´e omezuj´ı r˚ ust chyb v pr˚ ubˇehu v´ypoˇctu. Pˇr´ıklad: Soustava line´arn´ıch rovnic s matic´ı bl´ızkou k singul´arn´ı (ˇspatnˇe podm´ınˇen´a matice). Necht’ je d´ana ´uloha x+α y = 1 α x+y = 0 Necht’ vstupem je hodnota α a v´ystupem hodnota x. Pak 1 x= 1 − α2
a
Cp =
|δx| |x| |δα| |α|
Pˇri α2 → 1 je ´uloha ˇspatnˇe podm´ınˇen´a.
10
≃
dx α dα x
=
2 α2 |1 − α2 |
6
Numerick´ a stabilita
U nestabiln´ı metody se relativnˇe mal´e chyby v jednotliv´ych kroc´ıch v´ypoˇctu postupnˇe akumuluj´ı tak, ˇze dojde ke katastrof´aln´ı ztr´atˇe pˇresnosti numerick´eho ˇreˇsen´ı ´ulohy. U stabiln´ıch metod roste chyba v´ysledku s poˇctem krok˚ u N nejv´yˇse line´arnˇe (v √ ide´aln´ı, ale vz´acn´e situaci, kdy je znam´enko chyby n´ahodn´e, chyba roste ∼ N ). U nestabiln´ıch metod roste chyba rychleji, napˇr. geometrickou ˇradou ∼ q N , kde |q| > 1. Nestabilita metody vznik´a v d˚ usledku akumulace bud’ zaokrouhlovac´ıch chyb nebo chyb metody. Typicky se objevuje v rekurzivn´ıch algoritmech. Stabilita metody m˚ uˇze z´aviset na velikosti pouˇzit´eho kroku h. Nestabilita v d˚ usledku chyb metody se ˇcasto objevuje pˇri numerick´em ˇreˇsen´ı poˇc´ateˇcn´ıho probl´emu pro obyˇcejn´e a parci´aln´ı diferenci´aln´ı rovnice. 6.1
Pˇ r´ıklady nestabiln´ıch algoritm˚ u
1. Nestabiln´ı rekurze - uk´aˇzeme si na ponˇekud umˇel´em pˇr´ıpadu poˇc´ıt´an´ı mocnin ˇc´ısla Φ zvan´eho ”Zlat´y ˇrez” √ 5−1 Φ≡ ≃ 0, 61803398 2 Lehce uk´aˇzete, ˇze mocniny Φn splˇnuj´ı jednoduch´y rekursn´ı vztah Φn+1 = Φn−1 − Φn Protoˇze zn´ame Φ0 = 1 a Φ1 = 0, 61803398 mohli bychom zkusit poˇc´ıtat mocniny odeˇc´ıt´an´ım, coˇz je obvykle rychlejˇs´ı neˇz n´asoben´ı.
11
Výsledky výpoctu Φn
Hodnota, Chyba
1
0.5
0
−0.5
−1 0
Mocnina Výsledek rekurze Relativní chyba 10
20 Mocnina n
30
40
Obr´azek ukazuje, ˇze uveden´y postup zcela nepouˇziteln´y, pˇri jednoduch´e pˇresnosti dostaneme viditeln´e chyby v´ysledky uˇz od n = 16, kdy Φn ≃ 5.10−4. Pro n = 20 dostanu poprv´e z´aporn´y v´ysledek rekurze, a tedy rekurze uˇz nijak neaproximuje hodnotu mocniny. Nejdˇr´ıve vzroste relativn´ı chyba (chyba mˇen´ı znam´enko), pak se objev´ı z´aporn´e hodnoty Φn a nakonec zaˇcne dokonce r˚ ust absolutn´ı hodnoty Φn . Nestabilita se projev´ı i ve dvojit´e pˇresnosti, zaokrouhlovac´ı chyba nar˚ ust´a ale z menˇs´ı hodnoty a tak by se 1. z´aporn´y v´ysledek rekurze objevil pro n = 40. Pˇr´ıˇcina nestability √ je v tom, ˇze uveden´a rekurzn´ı formule m´a jeˇstˇe druh´e ˇreˇsen´ı Φ2 = −( 5 − 1)/2 < −1 < −Φ. Protoˇze rekurzivn´ı relace je line´arn´ı, absolutn´ı velikost zaokrouhlovac´ı chyby bude nar˚ ustat geometrickou ˇradou s kvocientem q = |Φ2| > 1. Protoˇze nav´ıc ˇreˇsen´ı kles´a, relativn´ı velikost zaokrouhlovac´ı chyby roste geometrickou ˇradou s kvocientem q ′ = |Φ2|/Φ > 1. Stejn´a rekurze by ale mohla b´yt pouˇzita pro v´ypoˇcet mocnin ˇc´ısla Φ2. Pro tento v´ypoˇcet je metoda stabiln´ı a pracovala by uspokojivˇe. 12
Uveden´y pˇr´ıklad byl umˇel´y, nicm´enˇe u mnoha speci´aln´ıch funkc´ı (napˇr. Besselovy funkce) se k v´ypoˇctu hodnoty funkc´ı pouˇz´ıvaj´ı podobn´e rekurzivn´ı relace, vˇzdy ovˇsem tak, aby metoda byla stabiln´ı. 2. Nestabiln´ı metoda pro v´ypoˇcet obyˇcejn´ych diferenci´aln´ıch rovnic Necht’ ˇreˇs´ıme obyˇcejnou diferenci´aln´ı rovnici 1. ˇr´adu y ′ = f (x, y) Na pˇr´ıkladu rovnice y ′ = −y s poˇc´ateˇcn´ı podm´ınkou y(0) = 1 ˇreˇsen´e ve smˇeru r˚ ustu promˇenn´e x uk´aˇzeme, ˇze dvoukrokov´a metoda 2. ˇr´adu y′ ≃
y(x + 2h) − y(x) = −y(x + h) 2h
je nestabiln´ı. Jde vlastnˇe o podobnou rekurzi √ jako v´yˇse a pro pomˇer q = y(x+h)/y(x) existuj´ı 2 ˇreˇsen´ı, q1 = −h+ 1 + h2 je v absolutn´ı hodnotˇe menˇs´ı neˇz 1 a odpov´ıd´a prvn´ım tˇrem ˇclen˚ um Taylorova rozvoje ˇreˇsen´ı √ 2 y = y(0) · exp(−x), druh´y koˇren q2 = −h − 1 + h je v absolutn´ı hodnotˇe vˇetˇs´ı neˇz 1 a zp˚ usobuje nestabilitu algoritmu. Na n´asleduj´ıc´ım grafu je porovn´ana celkov´a relativn´ı chyba uveden´e nestabiln´ı metody s chybou Eulerovy metody y(x + h) = y(x) + h · y ′ (x) = y(x) − h · y(x) Eulerova metoda se obvykle nepouˇz´ıv´a, nebot’ jde o metodu 1. ˇr´adu s velkou chybou metody, nicm´enˇe je pro uveden´y pˇr´ıpad stabiln´ı.
13
Diferenciální rovnice y´= − y (y(0)=1)
Hodnota, Chyba
1
0.5
0
−0.5
−1 0
exp(−x) Rel. chyba − nestabilní metoda Rel. chyba − Eulerova metoda 1
2
3
4
x Obr´azek ukazuje, ˇze na zaˇc´atku ˇreˇsen´ı je nestabiln´ı metoda vzhledem k relativnˇe mal´e chybˇe metody pˇresnˇejˇs´ı, ale postupn´y r˚ ust chyby pˇrivede nakonec ke katastrof´aln´ım chyb´am. Katastrof´aln´ım chyb´am nelze zabr´anit zkracov´an´ım kroku, uˇzit´ı dvojit´e pˇresnosti katastrofu pouze odd´al´ı. U stabiln´ı metody roste chyba s d´elkou intervalu nejv´yˇse line´arnˇe a chybu lze zmenˇsit zkracov´an´ım kroku. 3. Nestabiln´ı spline Pˇri interpolaci dat pomoc´ı kubick´eho splinu (lok´aln´ı interpolace kubick´ym polynomem se spojitou derivac´ı) je tˇreba zadat 2 podm´ınky (napˇr. hodnotu derivace funkce) v obou krajn´ıch bodech. Nespr´avnou a nestabiln´ı metodu dostaneme, pokud obˇe podm´ınky zad´ame v 1 z okrajov´ych bod˚ u. Pokud jako 2. podm´ınku v prvn´ım okrajov´em bodu zad´ame napˇr. jako 2. derivaci rovnou hodnotˇe druh´e derivace, kter´a vyˇsla pˇri stabiln´ım postupu, tj. zadan´ych 1. derivac´ıch v obou okrajov´ych bodech, obˇe u ´lohy jsou z ma14
5
tematick´eho hlediska zcela ekvivalentn´ı a v pˇr´ıpadˇe poˇc´ıt´an´ı s pˇresn´ymi ˇc´ısly bych dostal totoˇzn´y v´ysledek. Pokud vˇsak numericky poˇc´ıt´am s koneˇcnou d´elkou ˇc´ısel, zaokrouhlovac´ı chyba vˇsak pˇri postupn´em poˇc´ıt´an´ı od 1 okraje nar˚ ust´a a ˇreˇsen´ı zaˇcne mezi zadan´ymi body silnˇe oscilovat.
Interpolace 50 hodnot sin(x) kubickým splinem 1
Hodnota
0.5
0
−0.5 Stabilní metoda Nestabilní metoda −1 0
1
2
3 x
4
5
6
Na grafu jsou uk´az´any oscilace chybnˇe poˇc´ıtan´eho (nestabiln´ıho) kubick´eho splinu. Je vidˇet jen nˇekolik prvn´ıch extr´em˚ u, ostatn´ı jsou v absolutn´ı hodnotˇe pˇr´ıliˇs velk´e (aˇz 1013). Pˇri poˇc´ıt´an´ı v dvojit´e pˇresnosti se viditeln´e oscilace objev´ı pro x > 4.
15
7
Volba metody (algoritmu) • Z´ akladn´ım poˇ zadavkem je moˇ znost vyˇ reˇ sen´ı u ´ lohy s dostateˇ cnou 1 ˇ pˇ resnost´ı. Casto je sledov´ana konvergence , coˇz znamen´a schopnost vyˇreˇsit u ´lohu s libovolnˇe vysokou pˇresnost´ı (omezen´e jen zaokrouhlovac´ı chybou) pˇri kroku h → 0 (pˇri poˇctu operac´ı N → ∞). • Pˇri v´ybˇeru metody hraje roli i sloˇzitost algoritmu (poˇcet operac´ı nutn´ych k z´ısk´an´ı v´ysledku se zadanou pˇresnost´ı) a pamˇet’ov´e n´aroky. • Je k dipozici spolehliv´a implementace pˇr´ısluˇsn´e metody?
Numerick´ e knihovny • Pro drtivou vˇetˇsinu ´uloh jsou k dispozici procedury ve standardn´ıch knihovn´ach. Pokud ´uloha nen´ı trivi´aln´ı, neprogramuji ji s´am! • Vˇetˇsina knihoven je ve FORTRANU • Profesion´aln´ı knihovny jsou drah´e (b´yvaj´ı k dispozici na velk´ych poˇc´ıtaˇc´ıch) - nejzn´amˇejˇs´ı NAG, IMSL - ale nˇeketr´e ˇc´asti jsou implementov´any v programov´ych bal´ıc´ıch jako MAPLE nebo ORIGIN • Pro uk´azky budeme pouˇz´ıvat knihovny NUMERICAL RECIPES (je pˇr´ılohou knihy) - FORTRAN, C, Pascal, C++, Python, Java a existuje interface k MATLABu • - volnˇe (byt’ ˇcasto s omezen´ımi) dostupn´y software - vyhled´av´an´ı na http://gams.nist.gov, mnoho softwaru je na serverech NETLIB (www.netlib.org), se zrcadly napˇr. na http://netlib.no v Bergenu, Norsko nebo na http://netlib.sandia.gov v SNL, USA.
1
konvergenci budeme pˇresnˇeji definovat pro konkr´etn´ı u ´ lohy
16
