Skripta pracovní verze

me 1. a r-t łÆdek a z skÆme (A;b)
(iii) pro i = 2;3;:::;n odeŁteme v matici (A;b) (ar1=a11)-nÆsobek 1.łÆdku
od i-tØho łÆdku a z skÆme matici:
(A0jb0) =
0
BB
BB
@
a011 a012 ::: a01n
0 a022 ::: a02n
... ...
0 a0n2 ::: a0nn





b01
b02
...
b0n
1
CC
CC
A
:
Proces opakujeme ve 2. sloupci s r = 2;3;:::;n atd. Postupn z skÆme
z matice A horn trojœheln kovou matici.
4.2 Stabilita, podm n nost a anal za chyb
V souladu s œvodn kapitolou pova ujeme algoritmus (v poŁetn metodu,
numerick postup) łe„en systØmu Ax = b za stabiln , jestli e pro vypoŁtenØ
łe„en bx je
(A+ A)bx = b+ b; (4.3)
kde A (tzv. chybovØ matice k matici A) a b (tzv. vektor chyb pravØ
strany systØmu) jsou dostateŁn malØ.
PoznÆmky
1. Malost v „e uveden ch veliŁin se chÆpe ve smyslu jejich absolutn ch,
nebo relativn ch chyb po komponentÆch a pł p. s vyu it m jejich norem.
2. Napł. aplikace GEM na matici
A =
10 10 1
12
!
vede k
A(1) = A;A(2) =
10 10 1
0 1010
!
:
płitom zm na (tzv. røstov faktor) mezi A(1) a A(2) je 121010: Proto płi
łe„en systØmu s matic A GEM nelze oŁekÆvat malou chybovou matici
E a tedy metoda nen stabiln (systØm Ax = (1=3) mÆ łe„en x1 = 1,
x2 = 1, ale A(3) = 1=3 bx1 = 0, bx2 = 1 s płesnost 10 3).
23
3. Obecn , GEM bez v b ru pivota je nestabiln , GEM s ŁÆsteŁn m a
nebo œpln m v b rem pivota pł p. GEM aplikovanÆ na systØm se spe-
ciÆln matic (symetrickÆ, diagonÆln dominantn , pÆsovÆ...) stabiln .
Pł klad 4.1 GEM łe„me systØm:
2x1 + 4x2 x3 = 5
x1 + x2 3x3 = 9
4x1 + x2 + 2x3 = 9
tj. (Ajb) = (A(1)jb(1)) =
0
B@ 2 4 11 1 3
4 1 2




5
9
9
1
CA:
Pak
A(2)jb(2)) =
0
B@ 2 4 10 1 5=2
0 7 4




5
13=2
19
1
CA
a
A(3)jb(3)) =
0
B@ 2 4 10 1 5=2
0 0 43=2




5
13=2
129=2
1
CA
Nyn łe„ me systØm A(3)x = b(3), tj.
2x1 + 4x2 x3 = 5
x2 5=2x3 = 13=2
43=2x3 = 129=2
e„ me od posledn rovnice (zp tn chod):
x3 = 3; x2 = 1; x1 = 1:
V ta 4.3 Nech» A je regulÆrn matice a pro n jakou jej normu je
k Ak< 1kA 1k: (4.4)
Pak łe„en ^x = x0+ x0 systØmu (4.4) aproximuje łe„en x0 systØmu (4.1)
s chybou
k x0k
kx0j <
k(A)
1 k(A)k AkkAk
k bk
kbk +
k Ak
kAk
!
24
PoznÆmky
1. Je-li b = 0 (pravÆ strana systØmu je dÆna płesn ), nebo A = 0 (matice
A systØmu je dÆna płesn ), stejn jako v obecnØm pł pad - ovliv uje
relativn chybu łe„en Ł slo podm n nost matice A, tj. k(A).
2. Pł kladem „patn podm n nØ matice je tzv. Hilbertova matice
H =
0
BB
BB
B@
1 1=2 1=3 ::: 1=n
1=2 ...
... ...
1
n
1
n+1 :::
1
2n 1
1
CC
CC
CA
Pro n = 10 je k(A) := 3;21013 (s nn1), aŁ jHj> 1,
e je „patn podm n nÆ.
25
Pł klad 4.2 k(E) = 1, tj. E je dobłe podm n na.
Pł klad 4.3 Je-li A matice soustavy z poznÆmky 4.2.4, k(A) = 1112 =
12321, tj. A je „patn podm n na.
Bu nyn A matice s poruchou A, b vektor s poruchou b a spolu se
systØmem (4.1) łe„me i systØm
(A+ A)(x+ x) = b+ b (4.5)
4.2.1 IteraŁn metody
Pł mØ metody łe„en (vŁetn GEM) nejsou pł li„ vhodnØ pro velkØ matice A
(n> 104). Situaci mohou usnadnit iteraŁn metody.
Princip
spoŁ vÆ ve vyjÆdłen pøvodn ho systØmu Ax = b v ekvivalentn m tvaru
x = Tx+g (4.6)
płiŁem x0 je łe„en m systØmu (4.1) prÆv tehdy, je-li łe„en m i systØmu
(4.6). Analogicky s œvahami o systØmu (4.1), je-li x0 łe„en m systØmu (4.1) s
regulÆrn matic A, je-li x0 = A 1b, je-li x0 łe„en m systØmu (4.6) a (E T)
regulÆrn , pak
x0 = (E T) 1g:
je-li x0 2 Rn libovolnÆ poŁÆteŁn aproximace, pak posloupnost
n
x(k)
o1
k=0urŁenÆ vztahem
x(k+1) = Tx(k) +g; k = 0;1;::: (4.7)
se naz vÆ iteraŁn posloupnost a T iteraŁn matice.
Z rovnice (4.7) plyne, e
x(k+1) = T(Tx(k 1) +g) +g = ::: =
= Tk+1x0 + (Tk +Tk 1 +:::+T +E)g:
Proto łekneme, e matice T je konvergentn , jestli e limk!+1Tk = .
Pł klad 4.4 T =
1=2 0
1=4 1=8
!
je konvergentn , proto e
Tk = ::: =
0
@ (
1
2)
k 0
k
2(k+1)
1
2
k
1
A
26
a
lim
k!1
Tk =
V ta 4.4 NÆsleduj c tvrzen jsou ekvivalentn :
(i) matice T je konvergentn
(ii) limk!1

Tk

= 0 pro n jakou normu k:k
(iii) (T) < 1 (kde je tzv. spektrÆln polom r T)
(iv) posloupnost
n
x(k)
o1
k=0 urŁenÆ iteraŁn m procesem (4.7) pro ka dou po-
ŁÆteŁn aproximaci x0 2Rn konverguje a limk!1x(k) = x0 je łe„en m
systØmu (4.1)
PoznÆmky
1. kriteria pro zastaven v poŁtu po dosa en płesnosti > 0kx(k+1) x(k)kkx(k)k <
nebo

r(k+1)

(kAk

xk+1

+kbk) kde r(k+1) = Ax(k+1) b je tzv.
reziduum po (k + 1)-n iteraci.
2. RøznØ iteraŁn metody z skÆme volbou matice T.
Jacobiova iteraŁn metoda
Zapi„me A = D L U; kde
D =
0
BB
@
a11 0
...
0 ann
1
CC
A
L =
0
BB
BB
B@
0 0
a21 ...
... ... ...
an1 an;n 1 0
1
CC
CC
CA
U =
0
BB
BB
B@
0 a21 0
... ... ...
... a
n 1;n
0 0
1
CC
CC
CA
27
D je diagonÆln matice, L je doln trojœheln kovÆ matice s nulami na
diagonÆle a U je horn trojœheln kovÆ matice s nulami na diagonÆle.
Rovnici Ax = b zap „eme ve tvaru (D L U)x = b a transformujeme ji
na rovnici
Dx = (L+U)x+b
Za płedpokladu, e aii 6= 0;i = 1;:::;n je matice D regulÆrn a z płed-
choz rovnice plyne maticov tvar tzv. Jacobiovy iteraŁn metody (JIM),
e x = Tyx+g, kde TJ = D 1(L+U) a g = D 1b. Płitom
TJ = (tij)ni;j=1; kde tij =f0 pro i=j aij
aii pro i6=j
Z V ty 4.4 plyne pro JIM:
V ta 4.5 Posloupnost fx(k)g1k=0 generovanÆ rovnic
x(k+1) = TJx(k) +g; k = 0;1;:::
konverguje pro ka dou poŁÆteŁn aproximaci x02Rn prÆv tehdy, kdy (TJ) <
1:
Pł klad 4.5 JIM łe„me systØm Ax = b, kde
A =
0
B@ 10 2 2 1 10 2
1 1 10
1
CA; b = (6;7;8)T
Pak
TJ =
0
B@ 0 0;2 0;20;1 0 0;2
0;1 0;1 0
1
CA; g = (0:6;0:7;0:8)T
(TJ) = 0;285 < 1; kTJk1 = 0;4:
Metoda tedy konverguje, płiŁem pro x0 = (0;0;0)T je
x1 =
0
B@ 0;60;7
0;8
1
CA x2 =
0
B@ 0;900;92
0;93
1
CA x3 =
0
B@ 0;9700;976
0;982
1
CA x4 =
0
B@ 0;99180;9934
0;9958
1
CA
Pro odhad chyby płitom plat

x xk

1 kTJk
k
1
1 kTk1

x1 x0

1;
28
co v na„em pł pad pro k = 4

x x4

1 0;4
4
0;6 0;8 0;034:
SkutreŁnÆ chyba je v„ak men„ , v „e uveden m zpøsobem je odhadnuta. płesnØ
łe„en je toti x0 = (1;1;1)T a

x0 x(4)

18:10 3.
29
Kapitola 5
Interpolace
C lem interpolace je nahradit funkci f pro v poŁty røzn ch charakteristik
funkce f vhodn j„ funkc dostateŁn bl zkou funkci f.
ZÆsadn m problØmem je v b r vlastnost funkce (z jakØho prostoru) a
kriteria urŁuj c ho bl zkost obou funkc :
- za jsou brÆny polynomy, goniometrickØ funkce, funkce spojitØ, hladkØ,...
- za kriterium bl zkosti rovnost funkŁn ch hodnot ve vybran ch bo-
dech, pł p. derivac , dosa en minima ve vhodn konstruovan ch funk-
cionÆlech: minparam ( c(f; ));:::
5.1 PolynomiÆln interpolace
Bu f funkce, kterÆ v n+ 1 røzn ch bodech xi(i = 0;1;:::;n) nab vÆ hod-
not f(xi)(i = 0;1;:::;n). Najd me polynom Pn = Pn(x) stupn nejv „e n
takov , e
Pn(xi) = f(xi) (i = 0;1;:::;n) (5.1)
Body xi(i = 0;1;:::;n) naz vÆme uzly a polynom Pn s vlastnost (5.1)
interpolaŁn polynom.
V ta 5.1 Pro (n+ 1) dan ch dvojic Ł sel
(xi;f(xi)); i = 0;1;:::;n xi6= xj pro i6= j
existuje prÆv jeden polynom Pn stupe n:
Pn(x) = a0xn +a1xn 1 +:::+an 1x+an
takov , e plat (5.1)
30
Z teorie polynomø płitom plyne, e
Pn(x) = l0(x)f0 +l1(x)f1 +:::+ln(x)fn; (5.2)
kde
li(x) = (x x0):::(x xi 1)(x xi+1):::(x xn)(x
i x0):::(xi xi 1)(xi xi+1):::(xi xn)
a
fi = f(xi) (i = 0;1;:::;n):
Tento polynom naz vÆme Lagrangeøv interpolaŁn polynom.
Pł klad 5.1 Sestrojte Lagrangeøv interpolaŁn polynom, je-li dÆno:
xi 0 1 2 5
fi 2 3 12 147
e„en : V tomto pł pad je n = 3; hledÆme tedy polynom
P3(x) = 2(x 1)(x 2)(x 5)( 1)( 1)( 5) + 3(x 0)(x 2)(x 5)(1 0)(1 2)(1 5) +
+12(x 0)(x 1)(x 5)(2 0)(2 1)(2 5) + 147(x 0)(x 1)(x 2)(5 0)(5 1)(5 2) =
= x3 +x2 x+ 2:
PoznÆmky
1. Podobn lze konstruovat interpolaŁn polynom k funkci f z dan ch
uzlø.
2. InterpolaŁn polynom Pn mø e b t zapsÆn i v jin ch tvarech (napł.
Newtonøv, ¨eby„evøv, :::) vhodn j„ ch napł. ke studiu chyb interpo-
lace.
3. Odhad chyby interpolace fce f polynomem Pn na intervalu [a;b] za
płedpokladu, e funkce f mÆ na tomto intervalu derivace a do łÆdu
(n+ 1) vŁetn :
jf(x) Pn(x)j Mn+1(n+ 1)! (b a)
n+1
22n+1 ;
kde

f(n+1)(x)

Mn+1 8x2[a;b].
31
5.2 Interpolace v ekvidistantn ch uzlech
PłedpoklÆdejme nyn , e uzly xi (i = 0;1;:::;n) jsou ekvidistantn , tj.
existuje h> 0, tzv. krok, tak, e
xi = x0 +ih; i = 0;1;:::;n:
Konstrukc Lagrangeova interpolaŁn ho polynomu Pn pak z skÆme
V ta 5.2 Nech» uzly xi (i = 0;1;:::;n) jsou ekvidistantn s krokem h. Pak
Newtonøv interpolaŁn polynom Pn lze zapsat ve tvaru
Pn(x0 +th) = f0 +
nX
j=1
jf0
j! t(t 1):::(t j + 1);
kde x = x0 +th t2R je novÆ prom nnÆ,
jf0 =
nX
j=1
( 1)j i
j
i
!
fi
je tzv. j-tÆ diference v bod x0 a fi = f(xi); i = 0;1;:::;n
PoznÆmka
Podobn lze vyjÆdłit Pn s pomoc diferenc v jinØm bod , napł. xn
5.3 Interpolace pomoc splajnø
Płedchoz postup konstrukce interpolaŁn ho polynomu Pn v celØm uva o-
vanØm intervalu nen v dy v hodn . N kdy je œŁeln j„ aproximovat funkci
f napł. polynomy ni „ ho stupn po ŁÆstech. V takovØm pł pad mluv me
o polynomiÆln ch splajnech1, jejich nejdøle it j„ m reprezentantem jsou
tzv. kubickØ splajny.
Bu f funkce de novanÆ na intervalu [a;b] a a = x0 0 je
limn!1P(jf fnj< ) = 1;
tak mø eme hledanou nÆhodnou veliŁinou f urŁit na zÆklad statistickØho
vyhodnocen nÆhodn ch pokusø.
7.1 V poŁet urŁitØho integrÆlu
A1 Typ Jf = R10 f(x)dx, płiŁem 0 f(x) 1; x2[0;1]
OznaŁme mno ina bodø [x;y] takov ch, e x2 [0;1] a 0 y
f(x)
A jev: nÆhodn zvolen bod [x;y] 2 [0;1] [0;1] padne do pod-
mno iny
Z de nice a vlastnost tzv. geometrickØ pravd podobnosti plyne, e
pravd podobnost P(A), e nastane jev A
P(A) = m {ra m {ra[0;1] [0;1] =
Rb
a f(x)dx
1 =
Z b
a
f(x)dx:
Z druhØ strany je P(A) := mn , kde n je poŁet nÆhodn ch pokusø a m
je poŁet t ch, kdy bod [x;y] 2 , płiŁem aproximace Ł sla P(A) je
38
t m lep„ , Ł m pou ijeme v t„ ho mno stv pokusø (to je døsledek tzv.
zÆkona velk ch Ł sel).
Postup:
Pro i = 1;:::;n
1. s pou it m generÆtoru nÆhodn ch Ł sel (tabulek nÆhodn ch Ł sel)
zvolme v i-tØm pokusu dvojici [x(i);y(i)] takov ch, e x(i) 2 [0;1]
a y(i)2[0;1]
2. zva me, zda
y(i) f(x(i)) (7.1)
pokud ano, płipoŁt me 1 k poŁtu œsp „n ch pokusø
3. oznaŁmem - poŁet œsp „n ch pokusø a urŁeme hodnotuR10 f(x)dx :=
m
n .
Pł klad 7.1 UrŁeme płibli nou hodnotu
Z 1
0
x2dx(= [x
3
3 ]
1
0 =
1
3)
metodou Monte Carlo.
i x(i) y(i) f(x(i))6= x(i)2 y(i)6= f(x(i))?
1 0,5 0,7 0,25 0
2 0,6 0,3 0,36 1
3 0,3 0,3 0,9 0
4 0,6 0,4 0,36 0
5 0,9 0,2 0,81 1
6 0,1 0,6 0,01 0
7 0,6 0,8 0,36 0
8 0,8 0,1 0,64 1
9 0,5 0,5 0,25 0
10 0,3 0,8 0,09 0
J = mn = 310 := 0;3
39
A2 Typ Jg = Rba g(x)dx, płiŁem g6 0 a je ohraniŁenÆ na intervalu [0;1]:
OznaŁme
D = inffg(x)jx2[0;1]g
H = supfg(x)jx2[0;1]g
Złejm D