NUMERICKÁ MATEMATIKA
JAN MALÝ
1.
Aproximace £ísel
K aproximaci £ísel se pouºívá mnoho metod (limitní p°echod, itera£ní metody, °et¥zové zlomky), nejznám¥j²í je metoda rozvoje v °adu. Známe Taylorovy °ady elementárních funkcí, nap°.
ex = 1 + x + ln(1 + x) = x −
∞ X x3 x2 xn + + ··· = , 2! 3! n! n=0
∞ X x2 x3 xn + + ··· = (−1)n+1 , 2 3 n n=1
∞ X x2n+1 x3 x5 1 1+x ln =x+ + + ··· = , 2 1−x 3 5 2n + 1 n=0
arctg x = x −
arcsin x = x + Vzorec pro
(1 + x)a
−1 ≤ x ≤ 1,
∞ X a n a(a − 1) 2 a(a − 1)(a − 2) 3 x ,, x + x + ··· = n 2! 3! n=0
−1 < x < 1
∞ X 1 x3 1 · 3 x5 1 · 3 · · · · · (2n − 1) x2n+1 + + ··· = . 2 3 2·4 5 2 · 4 · · · · · (2n) 2n + 1 n=0
binomická °ada, dosazením a =
1 2 dostaneme vzorec pro usnad¬uje výpo£et aproxima£ního polynomu. Je to vzorec
se nazývá
Hornerovo schéma
−1 < x < 1,
∞ X x3 x5 x2n+1 (−1)n + + ··· = , 3 5 2n + 1 n=0
(1 + x)a = 1 + ax +
−1 < x ≤ 1,
√
1 + x.
a1 a2 a0 + a1 x + a2 x2 + · · · = a0 1 + x(1 + x . . . . a0 a1 Poznamenejme, ºe konvergence Taylorovy °ady na hranici konvergen£ního intervalu, pokud v·bec nastává, m·ºe být velmi pomalá. Pro efektivn¥j²í výpo£ty je t°eba vyuºít rozvoj v bod¥ blízkém nule (týká se vý²e uvedených vzorc·, coº jsou vesm¥s rozvoje se st°edem v nule).
Úlohy.
P°i srovnávání rychlosti aproximace je uºite£né p°edepsat chybu, zastavit výpo£et p°i dosaºené
p°esnosti a porovnat po£et provedených úkon· (nap°. se£tených s£ítanc· °ady). K odhadu chyby je moºné vyuºít známé hodnoty po£ítané konstanty, ale je to nesportovní. Kdyº v²ak zastavíme výpo£et podle vzdálenosti po sob¥ jdoucích aproximací, skute£ná chyba m·ºe být daleko v¥t²í.
1.
Pm 1 1 n k=0 m! . Prove¤te n ) £i na základ¥ vzorce srovnání, berte v úvahu pracnost výpo£tu (mocnina je opakované násobení) a p°esnost aproximace. Spo£t¥te aproximace £ísla
e
na základ¥ vzorce
(1 +
Pro výpo£et °ady pouºijte Hornerovo schéma.
2.
ln(1 + x) hodnotu x = a − 1 pro a > 2 nejde v·bec a pro a = 2, 3, . . . je pouºitelné ln a = − ln a1 = − ln(1 + x), kde x = 1−a a . 1+x a+1 Ale mnohem lep²í je ln a = ln pro x = . Zkuste pro a = 2, 3. 1−x a−1 3. K výpo£tu π lze pouºít rozvoje arctg x v Taylorovu °adu. Dosazení x = 1 dává jednoduchý vzorec, Jak aproximovat
a=2
ln a?
Dosadit do vzorce
je tragicky pomalé. Pro
ale konvergence je tragická. Lep²í je
π 1 1 = arctg 1 = arctg + arctg . 4 2 3 V literatu°e (nap°. [DI, kap. XII.5]) lze nalézt ra novanou formuli
π 1 1 = arctg 1 = 4 arctg − arctg . 4 5 239 1
Zkuste r·zné vzorce a porovnejte rychlost konvergence.
4.
Pomocí rozvoje funkce
π. 5.
arcsin x
v Taylorovu °adu v bodech
x=
1 2,
x=1
po£ítejte aproximaci £ísla
Porovnejte rychlost aproximace.
Po£ítejte
√
Porovnejte vzorec
q
1 a
√ a pomocí itera£ní metody yn+1 = 12 yan + yn a pomocí Taylorova rozvoje funkce 1 + x. 1 rychlost aproximace. (Je-li a ≥ 2, volíme x < 0 tak, aby bylo 1 + x = a a pouºijeme =
√1 ). a 2.
Na²im úkolem je nahradit funkci
f
Aproxima£ní polynomy
funkcí
Pf,
která je blízká
f
a snáze se po£ítá. Pomocí základních
po£etních operací lze p°esn¥ vy£íslit pouze racionální funkce. Teorie racionální aproximace p°esahuje rámec kurzu, zde se budeme zabývat polynomiální aproximací. P°edpokládáme znalost Taylorova polynomu:
Tn f (x) =
n X f (k) (x0 )
k!
k=0 kde st°ed
(x − x0 )k ,
x0 volíme nej£ast¥ji roven nule, ale nap°. níºe u odvození Newtonovy metody budeme pot°ebovat
obecný tvar. Aby Taylor·v polynom dob°e aproximoval zadanou funkci, je t°eba hlídat derivace vysokých °ád·. Dal²í nevýhodou je ubývající p°esnost s rostoucí vzdáleností od st°edu. P°esto nap°. pro aproximaci elementárních funkcí má zásadní význam a poskytuje nejelegantn¥j²í vzorce. V dal²ím budeme °e²it úlohu, jak sestrojit aproximující polynom, známe-li hodnoty funkce v kone£n¥ mnoha tzv. uzlových bodech. Uvaºujme d¥lení
a = x0 < x1 < · · · < xm = b intervalu
ha, bi
(body
xi
budeme nazývat uzlovými). Obecný vzorec bude
P f (x) =
(1)
m X
f (xj )pj (x),
j=0 kde polynomy 2.1.
p
a
pj
nezávisí na
f.
Lagrangeovy interpola£ní polynomy.
f
P°irozená my²lenka je volit polynomy
v uzlových bodech splývaly. Zvolme pro pevné
symbol
δij
i ∈ {0, 1, . . . , m}
funkci
fi
pi tak, aby hodnoty fi (xj ) = δij , kde
tak, ºe
(tzv. Kroneckerovo delta) je de nován p°edpisem
( δij =
0, i 6= j, 1, i = j,
Pak musí platit
δki = fk (xi ) = P fk (xi ) =
m X
fk (xj )pj (xi ) =
j=0 Hledáme tedy polynomy
pi
stupe¬ byl co nejniº²í. Máme
pi .
Jelikoº
pi
m X
δkj pj (xi ) = pk (xi ),
k, i = 0, 1, . . . , m.
j=0
pi (xj ) = δij , i, j = 0, 1, . . . , m. Volme je tak, aby m bod· mají být ko°eny hledaného polynomu m. Snadno se p°esv¥d£íme, ºe polynom m-tého stupn¥
tak, abychom dostali
m+1
d¥licích bod·, z toho
má být nenulový, stupe¬ je aspo¬
s p°edepsanými vlastnostmi je pro kaºdé
i
práv¥ jeden, a to
(x − x1 ) . . . (x − xm ) , (x0 − x1 ) . . . (x0 − xm ) (x − x0 )(x − x2 ) . . . (x − xm ) p1 (x) = , (x1 − x0 )(x1 − x2 ) . . . (x1 − xm ) ...
p0 (x) =
pm (x) =
(x − x0 ) . . . (x − xm−1 ) , (xm − x0 ) . . . (xm − xm−1 )
Lagrangeovy polynomy aproximují danou funkci p°esn¥ v uzlových bodech, ale mimo uzlové body jsou t¥ºko ukontrolovatelné. I pro omezenou funkci m·ºe maximální chyba jít k nekone£nu, jak ukazuje následující p°íklad. (Daná funkce je sice nespojitá, ale to je jen pro jednoduchost výpo£tu. Efekt lze pozorovat i na nekone£n¥krát diferencovatelných funkcích).
2
P°íklad 1.
f , která je na intervalu h−1, 1i 0 krom¥ po£átku, kde má hodnotu 1. Pro uzlové body nk , k ∈ {−n, . . . , −2, −1, 0, 1, 2, . . . , n}, má aproxima£ní polynom, který si pro te¤ ozna£me L2n f , tvar nx nx nx nx nx nx ... 1 + 1+ 1− 1− ... 1 − . L2n f (x) = 1 + n 2 1 1 2 n n−1 Uvaºujme n ≥ 3 a zvolme x > n . Potom nx > n − 1 a proto nx nx nx ... 1 + 1+ > 1 · 2 . . . 2 = 2n−1 . 1+ n 2 1 Pokusíme se o odhad maximální chyby p°i aproximaci funkce
v²ude rovna
Dále máme odhad
(nx − 1)(nx − 2) . . . (nx − n + 2)(nx − n + 1)(nx − n) nx nx nx 1− ... 1 − =− 1− 1 2 n n! (n − 2)(n − 3) . . . 1 · (nx − n + 1)(n − nx) ≥ n! (nx − n + 1)(n − nx) = . n(n − 1) Funkce
(nx − n + 1)(n − nx)
nabývá na intervalu
dostaneme
(ntn − n + 1)(n − ntn ) = odtud
L2n f (tn ) ≥ 2.2.
Bernsteinovy polynomy.
polynomy
pi
P (1 − f ) ≥ 0, pi ,
tn :=
2n 2n+1 , dosazením
n(n + 1) n(n − 1) > , 2 (2n + 1) (2n + 1)2
2n−1 →∞ (2n + 1)2
Nep°íjemné vlastnosti Lagrangeových polynom· jsou zp·sobeny tím, ºe
s maximální chybou v¥t²í neº
Polynomy
maxima v bod¥
hodn¥ kmitají. Výhodné je, kdyº aproxima£ní p°edpis (operátor
funkci také nezápornou funkci. Je-li navíc
Pf
h n−1 n , 1i
2.
P 1 = 1,
P
nem·ºe se stát, ºe by funkci
z (1)) p°i°adí nezáporné
f , −1 ≤ f ≤ 1,
p°i°adil
Jelikoº operátor je lineární, máme
P f = P 1 − P (1 − f ) ≤ P 1 = 1,
podobn¥
P f ≥ P (−1) = −1.
které aproximují nezáporná vstupní data, by m¥ly být téº nezáporné. Toho se dá docílit
h0, 1i a hledejme n-tého stupn¥, který má nejvý²e k ko°en· v intervalu (−∞, 0i, nejvý²e (n−k) ko°en· v intervalu h1, ∞), n¥kde nabývá hodnoty 1 a v ostatních bodech intervalu h0, 1i je co nejmen²í. Heuristickými
jedin¥ tím, ºe v²echny ko°eny leºí vn¥ daného intervalu. Omezme se nyní na interval polynom
úvahami dojdeme k záv¥ru, ºe povolený po£et reálných ko°en· je t°eba vyuºít a ºe tyto ko°eny je t°eba strkat co nejblíºe k bod·m
0, 1.
Poloºme tedy
qk (x) = ck xk (1 − x)n−k . Chytrá volba koe cient·
ck
vede k
Bernsteinovým polynom·m, které jsou dány vzorcem Bn f (x) =
n X
f
k n
bn,k (x),
k=0 kde
Pro
n k bn,k (x) = x (1 − x)n−k . k f (x) = 1
dostaneme z binomické v¥ty
Bn 1(x) =
n X n
k
k=0
n xk (1 − x)n−k = x + (1 − x) = 1.
Bernsteinovy polynomy sice neaproximují p°esn¥ v uzlových bodech, ale jinak dávají mnohem lep²í aproximaci neº Lagrangeovy polynomy. Pro danou funkci
Bn f (x)|
f
a maximální chybu
σn = maxx |f (x) −
lze dokázat, ºe
f
je spojitá
=⇒ σn → 0.
Tato vlastnost je teoreticky velmi významná a Bernsteinovy polynomy poskytují jeden z nejjednodu²²ích d·kaz· moºnosti takové tzv. stejnom¥rné aproximace. Chceme-li v²ak minimalizovat maximální chybu (p°i daném stupni polynomu), existují mnohem lep²í aproximace, dobré výsledky dávají tzv. eby²evovy
3
polynomy. Problém nejlep²í stejnom¥rné polynomiální aproximace je obtíºný a nedá se °e²it jednoduchou formulkou. 2.3.
L2 aproximace (Metoda nejmen²ích £tverc· I).
Místo minimalizace maximální chyby m·ºeme
chtít minimalizovat pr·m¥rnou chybu. Nejp°irozen¥j²í je minimalizovat pr·m¥r kvadrátu chyby (tzv. metoda nejmen²ích £tverc·). Ve skute£nosti nejde o metodu, ale o zadání úlohy, proto vhodn¥j²í by byl pojem
L2
aproximace. (íslo
b
Z
(f (x) − g(x))2 dx
1/2
a
L2 -vzdálenost funkcí f a g .) Hledáme nejmen²í hodnotu Z b (c0 , . . . , cn ) 7→ (f (x) − c0 − c1 x − c2 x2 − . . . − cn xn )2 dx.
se ve vy²²í matematice nazývá
funkce
a Tato funkce je kvadratický polynom
n
prom¥nných a minima nabývá tam, kde jsou nulové v²echny
parciální derivace. Derivováním dostaneme podmínky
Z
b
(f (x) − c0 − c1 x − c2 x2 − . . . − cn xn ) dx,
0= a
Z
b
x(f (x) − c0 − c1 x − c2 x2 − . . . − cn xn ) dx,
0= a
Z
b
x2 (f (x) − c0 − c1 x − c2 x2 − . . . − cn xn ) dx,
0= a
... Z 0=
b
xn (f (x) − c0 − c1 x − c2 x2 − . . . − cn xn ) dx,
a které se dají upravit na soustavu integrály 2.4.
Rb a
xk f (x) dx,
n lineárních rovnic o n neznámých c0 , . . . , cn . P°itom je pot°eba vy£íslit
coº m·ºe p·sobit potíºe.
Metoda nejmen²ích £tverc· II.
Namísto integrální vzdálenosti m·ºeme v metod¥ nejmen²ích
£tverc· m¥°it blízkost funkcí ve vybraných uzlových bodech. Máme-li d¥lení
ha, bi a £ísla κ0 , . . . , κm > 0 (váhy), m·ºeme hledat p = c0 + c1 x1 + · · · + cn xn stupn¥ nejvý²e n tak, aby výraz m X (2) κi (p(xi ) − f (xi ))2
intervalu
a = x0 < x1 < · · · < xm = b f : ha, bi → R polynom
k dané funkci
i=0 byl co nejmen²í. Podobn¥ jako v p°edchozím p°ípad¥, úloha op¥t vede na soustavu lineárních algebraic-
(n+1) neznámých c0 , . . . , cn+1 , kterou odvodíme derivováním (2) podle t¥chto prom¥nných. n ≥ m najdeme polynom, který má v bodech xi hodnoty f (xi ), pro n = m je tento polynom jednozna£n¥ ur£en a je to Lagrange·v interpola£ní polynom. Úloha za£ne být nová a zajímavá pro n < m. T°eba v p°ípad¥ n = 1 pro κi = 1 (lineární regrese) dostaneme rovnice X X nc0 + x i c1 = f (xi ) , kých rovnic o V p°ípad¥
i
X i 2.5.
Spline approximation.
xi c0 +
X
i
x2i
X c1 = xi f (xi ) .
i
i
V praxi nás nikdo nenutí, abychom pouºili jednotný aproxima£ní polynom
na celém daném intervalu. Tzv. aproximace pomocí splin· je zaloºena na my²lence, ºe aproximující polynom se m¥ní v kaºdém uzlovém bod¥. Zadáme-li hodnoty v sousedních bodech, m·ºeme je spojit lineárním polynomem. K dosaºení vy²²í p°esnosti je n¥kdy vhodné zadat tam i derivace jednoho nebo více °ád·. ím více derivací p°edepí²eme, tím pot°ebujeme vy²²í stupe¬ aproximujícího polynomu. Hledání aproxima£ního polynomu vede na soustavu lineárních algebraických rovnic. Uvaºujme funkci
f : ha, bi → R a d¥lení a = x0 < x1 < · · · < xm = b intervalu ha, bi. Lineární spline g a splývá s funkcí f v uzlových bodech. Tedy na intervalu hxi−1 , xi i je x − xi−1 xi − x f (xi−1 ) + f (xi ). g(x) = xi − xi−1 xi − xi−1
je po £ástech lineární funkce
4
g m¥la s funkcí f v uzlových bodech spole£né téº derivace do °ádu k , hledáme g 2k + 1 (coº vede na spliny lichého stupn¥). Kubický spline je na intervalu hxi−1 , xi i p(x) = c0 + c1 x + c2 x2 + c3 x3 t°etího stupn¥, jehoº koe cienty najdeme jako °e²ení soustavy
Chceme-li, aby funkce jako polynom stupn¥ polynom
£ty° lineárních rovnic
p(xi−1 ) = f (xi−1 ), p(xi ) = f (xi ), 0
p (xi−1 ) = f 0 (xi−1 ), p0 (xi ) = f 0 (xi ). Úlohy. 1.
Napi²te program na graf aproximujícího polynomu, na n¥mº by k zadaným uzlovým bod·m a hodnotám byl vid¥t sou£asn¥ Bernstein·v polynom a Lagrange·v interpola£ní polynom. Zkoumejte citlivost na zm¥nu zadané hodnoty v jednom bod¥.
2.
Napi²te program na graf aproximujícího polynomu, na n¥mº by bylo vid¥t srovnání Taylorova polynomu a polynomu nejlep²í
L2 -aproximace pro funkci ex . V²imn¥te si, jak chyba Taylorova polynomu L2 aproximace je rovnom¥rn¥ji
vzr·stá s rostoucí vzdáleností u po£átku, zatímco u polynomu nejlep²í rozprost°ena.
3.
Napi²te program na graf aproximujícího polynomu, který s danou funkcí na intervalu má spole£né hodnoty a derivace v krajních bodech. Porovnejte p°esnost s lineární interpolací nap°. pro
cos x
na
f (x) =
h0, 1i. 3.
Soustavy algebraických rovnic
Pro dané obecn¥ nelineární funkce
fi
prom¥nných
x1 , . . . , x n
máme °e²it soustavu rovnic
f1 (x1 , . . . , xn ) = 0, f2 (x1 , . . . , xn ) = 0, ... fn (x1 , . . . , xn ) = 0. Problém zapí²eme ve vektorovém tvaru
f (x) = 0.
(3) 3.1.
Metoda postupných aproximací.
Máme dánu rovnici (ve vektorovém tvaru)
x = g(x).
(4)
Na tento tvar lze p°evést rovnici (3) nap°. trikem
x = x + f (x).
Metoda postupných aproximací
je itera£ní metoda podle vzorce
xn+1 = g(xn ), V¥ta 2
(Banachova o pevném bodu).
Pokud existuje κ ∈ (0, 1) tak, ºe
Nech´ F ⊂ Rn je uzav°ená mnoºina a g : F → F je zobrazení.
|g(y) − g(x)| ≤ κ|y − x|,
(5)
potom °e²ení rovnice
(4)
x, y ∈ F,
existuje, je jednozna£né a metoda postupných aproximací k n¥mu konverguje.
Uv¥domme si, ºe p°edpoklady v¥ty nejsou jen formální a je t°eba ov¥°it, kdy jsou spln¥ny! Podívejme se nyní na odhad chyby. Máme
|x2 − x1 | = |g(x1 ) − g(x0 )| ≤ κ|x1 − x0 |, |x3 − x2 | = |g(x2 ) − g(x1 )| ≤ κ|x2 − x1 | ≤ κ2 |x1 − x0 |, ... |xn+1 − xn | = |g(xn ) − g(xn−1 )| ≤ κn |x1 − x0 |. Odtud pro
m>n |xm − xn | ≤ |xn+1 − xn | + · · · + |xm − xm−1 | ≤ (κn + · · · + κm−1 )|x1 − x0 |. 5
Pravou stranou m·ºeme odhadnout zbytkem geometrické °ady, takºe
|xm − xn | ≤ Limitní p°echod
m→∞
κn |x1 − x0 |. 1−κ
dává
|x − xn | ≤ Cκn , Tedy £ím
κ
kde
C=
|x1 − x0 | . 1−κ
je men²í, tím je rychlost konvergence v¥t²í.
Je dobré si uv¥domit, ºe metoda pracuje pro soustavy, ale demonstrovat si ji budeme pro jednoduchost
f (x) = 0 a p°edpokládejme, ºe f je diferencoI s nenulovou derivací. Poloºme g(x) = x + λf (x), kde λ 6= 0 si m·ºeme zvolit, a´ f (x) = 0 ⇐⇒ g(x) = x. Abychom dostali podmínku (5), pot°ebujeme a sta£í nám
stejn¥ na jedné rovnici o jedné neznámé. M¥jme rovnici vatelná na intervalu ho zvolíme jakkoli,
|g 0 | ≤ κ < 1.
Jelikoº
nejlep²í výsledky dostáváme pro
λ
blízké
g 0 (x) = 1 + λf 0 (x), 1 hodnotám − 0 f (x) .
Úlohy. 1.
x = g(x) metodou postupných aproximací s instruktivním√gra cg(x) = cos x demonstrujte konvergenci. Jak to dopadne pro g(x) = 1 − x,
Napi²te program na °e²ení rovnice kým výstupem. Na funkci
g(x) = 3x(1 − x)? 4.
Algebraické rovnice o jedné neznámé
Metoda p·lení interval·.
Nech´ f je spojitá funkce na intervalu ha, bi. e²me rovnici f (x) = 0. f (a) < 0 < f (b), m·ºeme postupovat takto: Chceme rekurentn¥ konstruovat body ak , bk tak, aby bylo f (ak ) ≤ 0 ≤ f (bk ), intervaly hak , bk i byly do sebe zano°ené a rychle se zmen²ovaly. Pak bude existovat spole£ná limita c = lim ak = lim bk , která bude ko°enem funkce f . Poloºme a0 = a, b0 = b. 1 Máme-li ak , bk , ozna£me ck = 2 (ak + bk ) (bod, který d¥lí interval hak , bk i nap·l). Jestliºe f (ck ) < 0, poloºme ak+1 = ck , bk+1 = bk . Jestliºe f (ck ) > 0, poloºme ak+1 = ak , bk+1 = ck . Jestliºe f (ck ) = 0, 4.1.
Zjistíme-li, ºe
m·ºeme si svobodn¥ vybrat z vý²e uvedených moºností, ale uv¥domme si, ºe ko°en uº jsme v tom p°ípad¥ na²li. Snadno spo£ítáme, ºe chybu metody lze odhadnout nerovností
|c − ak | ≤ bk − ak = 2−k (b − a), podobn¥ pro
|c − bk |.
Lze vymyslet rychlej²í metody, ale tato je naprosto spolehlivá.
4.2.
Newtonova metoda.
nici
f (x) = 0.
Nech´
f
je spojit¥ diferencovatelná funkce na intervalu
Zvolíme po£áte£ní aproximaci
polynomem prvého stupn¥ se st°edem v
x0
x0
hledaného ko°ene. Funkci
f
ha, bi.
e²me rov-
nahradíme Taylorovým
a dostaneme p°ibliºný vzorec
f (x) ∼ f1 (x) := f (x0 ) + f 0 (x0 )(x − x0 ). Funkci
f
nahradíme tímto jejím p°iblíºením
f1
a snadno °e²íme rovnici
f1 (x) = 0.
e²ení nazveme
x1
a
opakujeme proces:
f2 (x) := f (x1 ) + f 0 (x1 )(x − x1 ),
x2
je °e²ení rovnice
f2 (x) = 0.
Dostaneme obecný vzorec
f (xn ) = −f 0 (xn )(xn+1 − xn )
(6) a itera£ní metodu
xn+1 = xn −
f (xn ) f 0 (xn )
Pokusme se odhadnout chybu. P°edpokládejme, ºe na intervalu
00
|f (x)| ≤ L.
I,
kde hledáme °e²ení, je
f 0 (x) ≥ 1/K ,
Potom máme
P°edpokládejme, ºe
xn+1
i
|f 0 (y) − f 0 (x)| ≤ L|y − x|, x, y ∈ I. xn leºí v I . Z v¥ty o st°ední hodnot¥ dostaneme ξn
mezi
xn
s vyuºitím (6)
f (xn+1 ) = f (xn ) + f 0 (ξ)(xn+1 − xn ) = (f 0 (ξn ) − f 0 (xn ))(xn+1 − xn ), takºe
|f (xn+1 )| ≤ L|ξ − xn | |xn+1 − xn | ≤ L|xn+1 − xn |2 . 6
a
xn+1
tak, ºe
Dále
|xn+1 − xn | =
|f (xn )| ≤ K|f (xn )|, |f 0 (xn )|
takºe
|f (xn+1 )| ≤ K 2 L|f (xn )|2 . Poloºme
2
yn = f (xn ), C = K L
a ob¥ strany rovnice vynásobme
2
C.
Dostaneme
2
C|yn+1 | ≤ C |yn | , neboli pro
zn = Cyn
je
|zn+1 | ≤ |zn |2 . Poda°í-li se nám posloupnost
{xn }
|z0 | < 1,
rychlost konvergence je fantastická. Musíme ov²em je²t¥ pohlídat, aby se
udrºela v
I.
Nech´
x ¯
je °e²ení rovnice. Z v¥ty o st°ední hodnot¥ najdeme
ηn
tak, ºe
|xn − x ¯| |yn | = |f (xn ) − f (¯ x)| = |f 0 (ηn )| |xn − x ¯| ≥ , K takºe vzorec pro chybu je
|xn − x ¯| ≤ K|yn | ≤ Newtonova metoda konverguje velmi rychle, pokud nule a poda°í se nám uhodnout
x0
f
dostate£n¥ blízko
K |zn |. C
je dvakrát spojit¥ diferencovatelná,
x ¯.
f0
není blízké
Jinak ov²em m·ºe velmi snadno divergovat.
Úlohy. 1.
Napi²te programy na °e²ení rovnice
f (x) = 0
metodou postupných aproximací a Newtonovou meto-
dou s instruktivním gra ckým výstupem. Na p°íkladech (nap°.
f (x) = x2 − c)
porovnejte rychlost
konvergence obou metod. 5.
Soustavy lineárních rovnic
V této kapitole budeme °e²it soustavy rovnic
a11 x1 + · · · + a1n xn = b1 , a21 x1 + · · · + a2n xn = b2 ,
(7)
... am1 x1 + · · · + ann xn = bm ,
ve vektorovém zápisu
Px = q. Vektory x a q chápeme jako svislé a ztotoºníme-li se s maticemi o jednom Px = q interpretovat jako maticové násobení.
sloupci, m·ºeme po£etní úkon 5.1.
P°ímé metody.
Soustavy (7) se dají °e²it p°esn¥. Jediným d·vodem pro p°ibliºné numerické me-
tody je fakt, ºe pro velké soustavy je p°ímý výpo£et £asto del²í, neº p°ibliºný výpo£et vedoucí k dostate£né p°esnosti. P°ipome¬me tzv. Gaussovu elimina£ní metodu: P°ipome¬me, ºe Napi²me si je²t¥ matici o zprava vektor
m
°ádcích a
(n + 1)
P se nazývá matice soustavy. P p°ipí²eme
sloupcích, která vznikne, kdyº vedle matice
q (tzv. roz²í°ená matice soustavy). S °ádky provádíme ekvivalentní úpravy, coº jsou úpravy
zachovávající °e²ení soustavy. Jmenovit¥, m·ºeme prohodit po°adí °ádk·, nahradit °ádek jeho nenulovým násobkem, nebo nahradit °ádek jeho sou£tem s lineární kombinací ostatních °ádk·. M·ºeme téº vynechat °ádek, na kterém jsou samé nuly, nebo p°idat °ádek, který má na poslední pozici volitelné £íslo. (Z hlediska posuzování ekvivalence soustav takový °ádek vnímáme jako podmínku levá strana
∈ R,
takºe
neobsahuje ºádné omezení na °e²ení soustavy.) V kaºdém kroku m¥níme koe cienty a pravou stranu a m¥níme tím význam £ísel nyní obecný
k -tý
krok,
aij , bi .
Popí²eme
k = 1, . . . , n.
Rozli²íme dva p°ípady.
ai,k 6= 0. Zvolíme takové i. Máme-li více moºností, z hlediska i tak, aby prvek ai,q m¥l co nejv¥t²í absolutní hodnotu. Pokud jsme zvolili i 6= k , prohodíme i-tý a k -tý °ádek, takºe nyní je ak,k 6= 0. Vyd¥líme k -tý °ádek £íslem ak,k . Nyní je ak,k = 1. Pro v²echna i = k+1, . . . , m upravíme i-tý °ádek tak, ºe od n¥j ode£teme aik násobek k -tého °ádku. Nyní máme ai,k = 0 pro i > k . Tím je k -tý krok uzav°en. B) Pokud ai,k = 0 pro v²echna i ≥ k , posuneme v²echny °ádky od k -tého po£ínaje o jedno místo dol· a na k -tou pozici vsuneme °ádek, který má k -tý prvek 1, poslední prvek volitelné βk a jinak samé nuly. Po provedení v²ech n krok· zhodnotíme výsledek. Smaºeme p°ípadné nulové °ádky. Pokud nám zbude n¥jaký °ádek, který má na posledním míst¥ nenulu a jinak samé nuly (jelikoº akk = 1 pro k ≤ n, A) Nech´ existuje
i ∈ {k, . . . , m}
tak, ºe
minimalizace zaokrouhlovací chyby je nejlépe volit
7
takový °ádek musí být aspo¬
0=1
(n + 1)-vý),
je daná soustava ekvivalentní soustav¥ obsahující rovnici
a tudíº nemá ºádné °e²ení. V opa£ném p°ípad¥ má výsledná matice
n
°ádk·, je trojúhelníková,
na diagonále má jedni£ky a pod diagonálou nuly. Odpovídající soustava (s novými koe cienty a novou pravou stranou) má tvar
x1 + a12 x2 + a13 x3 + · · · + a1n xn = b1 , x2 + a23 x3 + · · · + a2n xn = b2 , ...
(8)
xn−1 + an−1,n xn = bn−1 , x n = bn , neboli
x n = bn , xn−1 = bn−1 − an−1,n xn , ...
(9)
x2 = b2 − a23 x3 − . . . − a2n xn , x1 = b1 − a12 x2 − a13 x3 − . . . − a1n xn . Nyní v tzv.
zp¥tném chodu
xn , . . . , x1 . V²imn¥me si, ºe na pravé stran¥ I je mnoºina index· k , pro varianty B . Potom v (8), (9) je bk = βk , tj. je to
z rovnic (9) postupn¥ vypo£ítáme
kaºdého °ádku se vyskytují jen ty neznámé, které uº jsou spo£ítané. Nech´ n¥º jsme v
k -tém
kroku algoritmu postupovali podle
volitelné £íslo. Pokud daná soustava má °e²ení, mnoºina v²ech °e²ení je a nní prostor, jehoº dimenze je po£et prvk· mnoºiny
I.
V tzv. regulárním p°ípad¥ je
I
prázdná mnoºina, do výpo£tu °e²ení nevstupují
volitelná £ísla a °e²ení je práv¥ jedno.
Itera£ní metody obecn¥.
5.2.
Matici soustavy si ozna£me
A,
n = m, tj. matice soustavy Ax = b upravíme na tvar
Budeme uvaºovat jen p°ípad
pravou stranu
b.
Rovnici
je £tvercová.
x = Px + q.
(10)
Nejprimitivn¥j²í zp·sob, jak to ud¥lat, je volba
P = I − A, kde
I
q = b,
zna£í jednotkovou matici, alle nikde není °e£eno, ºe máme pouºít zrovna tento p°edpis. Úlohu pak
°e²íme postupnými aproximacemi
xn+1 = Pxn + q,
(11)
n = 0, 1, . . .
Tato metoda nemusí konvergovat. Na pomezí lineární algebry a analýzy stojí d·leºitá v¥ta, která je maticovou analogií v¥ty o konvergenci geometrické °ady. Nejprve ale musíme p°ipomenout n¥které pojmy spektrální teorie. ekneme, ºe £íslo I kdyº
P
λ je vlastní £íslo
matice
P, jestliºe matice P − λI je singulární (tj. není invertibilní). λ ∈ C.
je reálná matice, je bohuºel t°eba se zaobírat i komplexními vlastními £ísly. Nech´
Následující podmínky jsou ekvivalentní: (i)
λ
P, Px = λx má nenulové °e²ení, n existuje vektor q ∈ R tak, ºe rovnice Px = λx + q hodnost matice P − λI je men²í neº n, det(P − λI) = 0. je vlastní £íslo
(ii) rovnice (iii) (iv) (v)
nemá °e²ení,
Následující v¥tu lze dokázat metodami algebry nebo metodami komplexní analýzy, oba d·kazy jsou t¥ºké.
P n Nech´ pro kaºdé vlastní £íslo λ matice P platí |λ| < 1. Potom °ada ∞ n=0 P je konvergentní P ∞ a sou£et n=0 Pn je inverzní matice k I − P. Itera£ní metoda (11) pak konverguje pro kaºdou pravou stranu q a kaºdou volbu po£áte£ní podmínky x0 .
V¥ta 3.
8
5.3.
Gauss-Seidelova metoda.
Abychom se vyhnuli index·m, itera£ní metodu budeme psát ve tvaru
y = Px + q, tedy místo
xn
budeme psát prost¥
M¥jme soustavu kde
L
Ax = b.
x a místo xn+1 . budeme psát y. A si m·ºeme rozepsat jako A = L + D + U (lower-diagonal-upper), D prvky na diagonále a U prvky nad diagonálou. Rovnici
Matici
jsou prvky pod diagonálou,
(L + D + U)x = b si upravíme na
(L + D)x = −Ux + b, x = −(D + L)
−1
neboli
Ux + (D + L)−1 b
Gaussova-Seidelova metoda je de nována pomocí formule
y = −(D + L)−1 Ux + (D + L)−1 b Zdálo by se, ºe kv·li pouºití Gauss-Seidelovy metody je t°eba invertovat. matici
D + L,
ale p°ekvapiv¥
ne. Metodu si p°epí²eme zp¥t do tvaru
(L + D)y + Ux = b, tedy odpovídající soustava je
a11 y1 + a12 x2 + a13 x3 + · · · + a1n xn = b1 , a21 y1 + a22 y2 + a23 x3 + · · · + a2n xn = b2 ,
(12)
... a21 y1 + a22 y2 + a23 y3 + · · · + ann yn = bm ,
e²íme-li úlohu postupn¥ shora, v kaºdém °ádku je tedy jen rovnice o jedné nové neznámá, která se snadno spo£ítá. Následující v¥ta se zabývá symetrickým p°ípadem
L = −T, U = −T∗ .
Minusy jsou tam proto, ºe
se s tím pak lépe po£ítá. Gauss-Seidelova metoda má celkem solidní pom¥r mezi pracností výpo£tu a konvergen£ními vlastnostmi.
Nech´ D je diagonální a T je dolní trojúhelníková s nulovou diagonálou. Je-li A = D − T − T∗ pozitivn¥ de nitní a P = −(D − T)−1 T∗ , pak pro v²echna vlastní £ísla λ matice P platí |λ| < 1. GaussSeidelova metoda pro Ax = b tedy v tomto p°ípad¥ konverguje. V¥ta 4.
D·kaz.
M¥jme
x ∈ Cn ,
ozna£me
y = −Px, u = x − y.
Máme
(D − T)y = T∗ x.
(13) k ob¥ma stranám (13) p°i£teme
−T∗ y
a dostaneme
Ay = T∗ x − T∗ y = T∗ u.
(14) K ob¥ma stranám (13) p°i£teme
Ax + (T − D)y
a dostaneme
Ax = (D − T)u.
(15) Rovnosti (14) a (15) vynásobíme
u,
kaºdou z jiné strany, a se£teme. Dostaneme
∗
u · Ax + Ay · u = u · T u + (D − T)u · u = Tu · u + Du · u − Tu · u = Du · u.
(16)
Úpravou levé strany s vyuºitím symetrie
A
dostaneme
u · Ax + Ay · u = x · Ax − y · Ax + Ay · x − Ay · y = Ax · x − Ay · y.
(17)
Z (16) a (17) máme
Ax · x − Ay · y = Du · u.
(18)
λ je vlastní £íslo P a x je odpovídající vlastní vektor. λ 6= −1, nebo´ jinak by bylo y = x, z rovnice (13) by plynulo
Nech´
Potom
y = −λx
a
u = (1 + λ)x.
Ax = (D − T)x − T∗ x = 0 a matice
A
by nemohla být pozitivn¥ de nitní. Podle (18) je
(1 − |λ|2 )Ax · x = Ax · x − λAx · λx = Ax · x − Ay · y = Du · u (19)
= (1 + λ)Dx · (1 + λ)x = |1 + λ|2 Dx · x. 9
Máme
Na pravé stran¥ (19) je kladný výraz, tedy na levé stran¥ také, takºe
Poznámka 5.
|λ| < 1.
Qx = d s regulární maticí soustavy. Rovnici si vynásobíme maticí Q∗ Q Qx = Q∗ d. Matice A = Q∗ Q je pak positivn¥ de nitní a m·ºeme pouºít Gauss-
M¥jme úlohu
zleva a dostaneme
∗
Seidelovu metodu. Zbývá jen singulární p°ípad, který nám tak snadno konvergentní metodu nenabízí. To je celkem pochopitelné, má-li úloha nekone£n¥ mnoho °e²ení (v lep²ím p°ípad¥!), pak p°íroda neví, které z nich si má vybrat.
Úlohy. 1. 2.
Napi²te program na °e²ení soustavy algebraických rovnic Gaussovou elimina£ní metodou. Napi²te program na °e²ení soustavy algebraických rovnic Gauss-Seidelovou metodou. 6.
Numerická integrace
Máme-li spo£ítat integrál
b
Z
f (x) dx, a v podstat¥ spole£ná my²lenka v²ech metod je nahradit funkci
f
blízkou funkcí a tuto zintegrovat. K tomu
m·ºeme pouºít výsledky kapitoly 2. Zde uvedeme metody, blízké my²lence de nice Riemannova integrálu. 6.1.
Obdélníkové metody.
intervalu
hxi−1 , xi i
a = x0 < x1 < · · · < xm = b intervalu ha, bi a v kaºdém ξi . V obdélníkové metod¥ nahradíme integrál sou£tem
M¥jme d¥lení
m¥jme je²t¥ zvolen bod
O=
m X
f (ξi )(xi − xi−1 ).
i=1
hxi−1 , xi i tedy funkci f nahradíme konstantou f (ξi ) a po£ítáme obsah obdélníka o rozm¥rech (xi − xi−1 ) × f (ξi ). Speciální volbou ξi dostaneme: Levou obdélníkovou metodu: ξi = xi−1 , pravou obdélníkovou metodu: ξi = xi , centrovanou obdélníkovou metodu: ξi = 21 (xi−1 + xi ). Na intervalu
Odhadn¥me nyní chybu pravé obdélníkové metody.
Nech´ funkce f na intervalu h0, 1i spl¬uje |f 0 | ≤ K . Potom, pouºijeme-li rovnom¥rné d¥lení na n interval·, chyba εn pravé obdélníkové metody je odhadnuta V¥ta 6.
K . 2n
εn ≤
(20)
D·kaz.
Na jednotlivém intervalu máme odhad
|f (x) − f (xi )| ≤ K(xi − x), takºe
Z
xi
Z f (x) dx − f (xi )(xi − xi−1 ) ≤
xi−1
xi
Z |f (x) − f (xi )| dx ≤
xi−1
xi
Z K(xi − x) dx =
xi−1
Se£tením p°es v²echny d¥licí intervaly dostaneme (20). 6.2.
Lichob¥ºníková metoda.
0
(−Kx) dx = 1 −n
K 2n2
a = x0 < x1 < · · · < xm = b intervalu ha, bi. hxi−1 , xi i funkci f nahradíme lineárním polynomem, který v krajjako f . a po£ítáme obsah lichob¥ºníka o vrcholech [xi−1 , 0], [xi , 0],
M¥jme op¥t d¥lení
V lichob¥ºníkové metod¥ na intervalu ních bodech nabývá stejné hodnoty
[xi , f (xi )], [xi−1 , f (xi−1 )].
Pouºijeme-li vzorec pro obsah lichob¥ºníku, dostaneme
L=
m X 1 i=1
2
[f (xi ) + f (xi−1 )](xi − xi−1 ).
Odhadn¥me nyní chybu metody.
Nech´ funkce f na intervalu h0, 1i spl¬uje |f 00 | ≤ K . Potom, pouºijeme-li rovnom¥rné d¥lení na n interval·, chyba εn lichob¥ºníkové metody je odhadnuta V¥ta 7.
(21)
εn ≤
K . 4n2
10
D·kaz.
Nech´
g
f . Funkce g ma v bodech xi−1 a xi f . Poloºme h = f − g , s = xi−1 , t = xi . Potom h(s) = h(t) = 0. Podle v¥ty o st°ední ξ ∈ (s, t) tak, ºe h0 (ξ) − 0. Pro derivaci funkce h tedy máme odhad
je po £ástech lineární funkce, kterou nahrazujeme
stejné hodnoty jako hodnot¥ existuje
|h0 (x)| ≤ K|x − ξ| ≤ Funkce
h
K . n
je tedy odhadnuta
K |x − c|, n
|h(x)| ≤
c = s, t.
a
Z
t
1 2n
Z |h(x)| dx ≤ 2
s
0
K K x dx ≤ 3 . n 4n
Se£tením p°es v²echny d¥licí intervaly dostaneme (21).
Poznámka 8.
Lichob¥ºníková metoda je tedy pro hladké funkce o °ád p°esn¥j²í neº pravá (£i levá)
metoda obdélníková. Je správné se na toto porovnání podívat kriticky. Co kdyº jsme ve v¥t¥ 6 pouºili
f (x) = x se snadno p°esv¥d£íme, ºe integrál n+1 1 1 a pravá obdélníková metoda dává aproximace , takºe chyba je °ádov¥ 2 2n n. Máme-li rovnom¥rné d¥lení nap°. intervalu h0, 1i na n interval·, snadnou úpravou dostaneme pro n-tý
ne²ikovný odhad, který se dá vylep²it? Na p°íkladu funkce je
sou£et tvar
1 1 f (0) + f n 2
Ln =
1 n
+f
2 n
+ ··· + f
n−1 n
+ 12 f (1)
zatímco obdélníková metoda je
On =
1 f n
1 n
2 n
+f
+ ··· + f
n−1 n
+ f (1) .
Pracnost výpo£tu je tedy prakticky stejná. Vidíme, ºe zdánliv¥ zanedbatelná úprava metody m·ºe p°inést zna£né zlep²ení p°esnosti. Poznamenejme je²t¥, ºe centrovaná obdélníková metoda dává srovnateln¥ dobrou p°esnost jako metoda lichob¥ºníková, otázka je, zda se nepo£ítají hodnoty funkce snáze v uzlových bodech neº uprost°ed mezi nimi. 6.3.
Simpsonova metoda.
intervalu
a = x0 < x1 < · · · < xm = b intervalu ha, bi. V kaºdém ξi = 21 (xi−1 + xi ). Potom integrál aproximujeme výrazem
M¥jme op¥t d¥lení
hxi−1 , xi i ozna£me ξi
jeho st°ed, tedy
m 1 X 2 1 (xi − xi−1 ) f (xi−1 ) + f (ξi ) + f (xi ) . 6 3 6 i=1 Výsledek odpovídá tomu, ºe na kaºdém intervalu nahradíme funkci souhlasí s funkcí
f
v bodech
xi−1 , xi
a
ξi
f
kvadratickým polynomem, který
a tuto po £ástech kvadratickou funkci zintegrujeme. Tato
metoda je je²t¥ o dva °ády p°esn¥j²í neº metoda lichob¥ºníková.
Úlohy. 1.
Napi²te program na numerickou integraci obdélníkovou a lichob¥ºníkovou metodou a porovnejte p°esnost dosaºených výsledk·.
7.
Numerické °e²ení oby£ejných diferenciálních rovnic
V této kapitole budeme vy²et°ovat diferenciální rovnici
y 0 = ϕ(x, y),
(22) kde
ϕ
je daná spojitá funkce prom¥nných
budeme rozum¥t diferencovatelnou funkci
x ∈ I , y ∈ R, I = ha, bi. e²ením rovnice f : I → R, která v kaºdém x ∈ I spl¬uje
(22) na intervalu
I
f 0 (x) = ϕ(x, f (x)). (v krajních bodech máme na mysli jednostranné derivace). V bod¥ podmínku (23)
f (a) = y0 . 11
a
si navíc zadáme tzv. po£áte£ní
7.1.
Eulerova metoda.
Uvaºujme d¥lení
a = x0 < x1 < · · · < xm = b intervalu
ha, bi.
P°ibliºné °e²ení rovnice (22), (23) budeme hledat tak, ºe derivaci nahradíme diferencí.
g : I → R, která bude lineární na intervalech d¥lení a uvnit° kaºdého intervalu g 0 (x) = ϕ(xi , g(xi )). Tedy dopustíme se té chyby, ºe ϕ x, ale v nejbliº²ím d¥licím bodu zleva xi . Ozna£me yi = g(xi ). Protoºe g je
Sestrojíme spojitou funkci
(xi , xi+1 )
bude spl¬ovat diferenciální rovnici
nebudeme vy£íslovat v lineární na
(xi , xi+1 ),
spl¬uje zde rovnici
g 0 (x) = Máme jiº zadanou po£áte£ní podmínku
yi+1 − yi g(xi+1 ) − g(xi )) = . xi+1 − xi xi+1 − xi y0 .
Známe-li
y0 , . . . , yi ,
z rovnice
yi+1 − yi = ϕ(xi , yi ) xi+1 − xi vypo£ítáme
yi+1 = yi + ϕ(xi , yi )(xi+1 − xi ). Odhad chyby budeme demonstrovat na rovnom¥rném d¥lení intervalu
h0, 1i
a autonomní rovnici
y0 =
ψ(y). V¥ta 9.
Nech´ funkce ψ : R → R na h0, 1i × R spl¬uje |ψ 0 | ≤ L.
|ψ| ≤ K,
Potom, pouºijeme-li rovnom¥rné d¥lení na n interval·, chyba εn Eulerovy metody (maximální odchylka |f (x) − g(x)|, kde f je skute£né °e²ení) je odhadnuta εn ≤
(24)
D·kaz.
KeL n
0 = x0 < · · · < x n = 1 , |g (x)| ≤ K , tedy
Uvaºujeme d¥lení
(xi , xi+1 )
máme
kde
xi =
0
i n . Ozna£me
zi = f (xi ).
Na intervalu
|g(x) − g(xi )| ≤ K(x − xi ), |f (x) − yi | ≤ |f (x) − g(x)| + |g(x) − g(xi )| ≤ |f (x) − g(x)| +
K , n
a
|f 0 (x) − g 0 (x)| = |ψ(f (x)) − ψ(yi )| ≤ L|f (x) − yi )|
≤ L|f (x) − g(x)| +
KL n
Ozna£me
K K h(x) = ln |f (x) − g(x)| + − ln n n Potom
|h0 | =
|f 0 (x) − g 0 (x)| ≤ L, |f (x) − g(x)| + K n
(to je velmi nep°esné, protoºe derivace ob£as nemusí existovat, ale dá se to spravit) a
h(0) = 0.
Tedy
|h(x)| ≤ Lx, neboli
|f (x) − g(x)| +
K n
K n
≤ eLx ≤ eL ,
takºe
|f (x) − g(x)| ≤
KeL n
12
Itera£ní metoda.
7.2. k
h.
Nejjednodu²²í diferenciální rovnice je
y 0 = h(x),
°e²ením je primitivní funkce
e²íme-li úlohu Eulerovou metodou, je to jako kdybychom provád¥li numerickou integraci levou
obdélníkovou metodou. Víme, ºe lep²í výsledky dává metoda lichob¥ºníková. Analogie lichob¥ºníkové metody na obecnou rovnici (25) Známe-li
yi ,
y 0 = ϕ(x, y) je metoda zaloºená na vzorci 1 yi+1 − yi ϕ(xi , yi ) + ϕ(xi+1 , yi+1 ) . = xi+1 − xi 2
t¥ºko m·ºeme z rovnice (25) p°ímo spo£ítat
yi+1 ,
protoºe rovnice o neznámé
yi+1
je impli-
citní. M·ºeme ale pouºít itera£ní numerickou metodu pro °e²ení algebraické rovnice a tím ur£it p°ibliºn¥
yi+1 , které je lep²í neº bychom dostali z Eulerovy metody. Najít optimální pom¥r mezi tím, kolik provád¥t iterací p°i hledání yi+1 a jak jemné d¥lení intervalu pouºít, je velké um¥ní. Tato úvaha je jen náznakem jak mohou za£ít snahy o ú£inn¥j²í metody °e²ení diferenciálních rovnic.
Úlohy. 1.
Napi²te program na numerické °e²ení diferenciální rovnice Eulerovou metodou s gra ckým výstupem.
13
