11
Interpolace, aproximace Metoda nejmensˇ´ıch cˇtvercu˚
11.1 Interpolace Mějme body [xi , yi ], i = 0, 1, . . . , n − 1. Cílem interpolace je najít funkci f (x), jejíž graf prochází všemi těmito body, tj. f (xi ) = yi . Poté lze pomocí nalezeného funkčního předpisu počítat hodnoty pro libovolné x, samozřejmě s jistou chybou. Nejčastější interpolační metody: Po částech lineární interpolace, kdy danými body proložíme lomenou čáru. Tj. dva sousední body jsou spojeny úsečkou, pro výpočet pak můžem užít oblíbenou trojčlenku. Polynomická interpolace, kdy danými body proložíme polynom. Dvěma body se dá vést přímka, tj. polynom stupně 1, třemi parabola, tj. polynom stupně 2 atd. Máme-li n bodů, můžeme jimi proložit polynom stupně n − 1: y = an−1 xn−1 + an−2 xn−2 + · · · + a1 x + a0 Dosadíme-li do tohoto polynomu [xi , yi ], i = 1, . . . , n, získáme soustavu n rovnic o n neznámých: ⎛ n−1 n−2 ⎞ ⎛ ⎞ ⎛ ⎞ x0 x0 . . . 1 an−1 y1 ⎜xn−1 xn−2 . . . 1⎟ ⎜an−2 ⎟ ⎜ y2 ⎟ 1 ⎜ 1 ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ .. .. .. ⎟ × ⎜ .. ⎟ = ⎜ .. ⎟ .. ⎝ . ⎠ ⎝ . . . . ⎠ ⎝.⎠ n−1 n−2 xn−1 xn−1 . . . 1 a0 yn Matice této soustavy se nazývá Vandermondova matice. Zásadní nevýhodou této metody je skutečnost, že budeme-li se snažit proložit více body (asi sedmi a více) jediný polynom, skončíme velmi pravděpodobně (pokud jsme data nezískali skutečně z nějakého polynomu) tak, že kvůli velmi výrazným a neopodstatněné zákmitům grafu bude polynom k ničemu (vizte obr. 1) Metoda je vhodná pouze při menším počtu bodů, kdy dává dobré výsledky.
1
Obr. 1: Polynom stupně 8 I při menším počtu bodů je výpočet koeficientů zatížen zaokrouhlovacími chybami (s Vandermondovou maticí se špatně počítá). Existují však i jiné postupy, jak polynom vypočítat. Newtonova interpolace: Newtonův interpolační polynom má následující tvar: Nn (x) = a0 +a1 (x−x0 )+a2(x−x0 )(x−x1 )+· · ·+an (x−x0 )(x−x1 ) . . . (x−xn−1 ) Koeficienty získáme postupným dosazováním argumentů, kterýžto postup lze algoritmizovat. Lagrangeova interpolace: Interpolační polynom L(x) =
n−1
yj j (x) = y0 0 (x) + y1 1 (x) + · · · + yn−1n−1 (x)
j=0
je lineární kombinací Lagrangeových bázových polynomů j (x) =
n−1 i=0, i=j
(x − xj−1 ) (x − xj+1 ) (x − xn−1 ) x − xi (x − x0 ) ··· ··· . = xj − xi (xj − x0 ) (xj − xj−1 ) (xj − xj+1 ) (xj − xn−1 )
Pro bázové polynomy j (x) platí: j (xi ) =
1 pro i = j 0 pro i = j
Existují i další metody, k jejich rozvoji vedla historicky snaha zjednodušit výpočetní náročnost a získat pohodlně interpolační polynom v případě doplnění dalších bodů k již použitým pro nalezení polynomu. Samozřejme, že všechny postupy vedou k témuž polynomu, ovšem při výpočtu různými metodami můžeme získat různé výpočty vlivem zaokrouhlovacích chyb. 2
Po částech polynomická interpolace, kdy např. třemi sousedními body proložíme parabolu. Interpolace splajnem [splajn, postaru často psáno spline], kdy mezi každými dvěma body je výsledná křivka popsána jiným funkčním předpisem fi (x). To umožňuje dosáhnout hladké funkce, tj. mající spojité derivace (nikde nejsou zuby). Často se užívá kubický splajn, kdy mezi každými dvěma body máme kubickou funkci fi (x) = ai + bi (x − xi ) + ci (x − xi )2 + di (x − xi )3 . Předpokládejme že máme n + 1 bodů (číslujeme 0, 1, . . . , n, tj. n intervalů), a tedy i n kubických polynomů, z nichž každý je popsán 4 koeficienty. Dohromady máme tedy 4n neznámých a musíme tedy mít 4n podmínek pro výpočet všech parametrů. Ty získáme takto: (1) jednostranné limity jednotlivých polynomů v bodech xi nám dají 2n rovnic, tedy: (1a) n podmínek typu lim fi (x) = yi x→xi +
pro i = 0, . . . , n − 1, kde fi (x) je odpovídající kubická funkce na i-tém intervalu (počítáno od nuly), (1b) n podmínek typu lim fi−1 (x) = yi x→xi −
pro i = 1, . . . , n, kde fi (x) je odpovídající kubická funkce na i-tém intervalu (počítáno od nuly), (2) z požadavku hladkosti nám vyplyne 2n − 2 podmínek typu (x) = lim fi (x) lim fi−1
x→xi −
x→xi +
(x) = lim fi (x) lim fi−1
x→xi −
x→xi +
pro i = 1, . . . , n − 1 (bez koncových krajních bodů), kde fi (x) je odpovídající kubická funkce na i-tém intervalu (počítáno od nuly). (3) pro dva koncové krajní body se obvykle bere podmínka tzv. přirozeného splajnu, tedy že do koncového bodu dobíhá splajn jako rovná čára, tedy lim f0 (x) = lim fn (x) = 0.
x→x0 +
x→xn −
Na pomoc si samozřejmě vezmeme kvalitní matematický softvér.
3
Obr. 2: Kubický splajn Mezi výhody kubického splajnu patří: (1) lze jednoduše přidávat další a další body tvořící složitější a složitější křivky při zachování přesnosti a proveditelnosti výpočtu (srovnejte si s polynomy vysokého stupně), (2) výsledná funkce je hladká, nikde žádné zlomy, (3) případný zákmit lze odstranit přidáním dalšího bodu.
11.2 Aproximace Mějme body [xi , yi ], i = 1, . . . , n. Argumenty xi jsou určeny deterministicky (tedy měřeny bezchybně) a yi jsou nezávislé náhodné veličiny (nabývají hodnot s jistou pravděpodobností) s normálním rozdělením se střední hodnotou 0 a rozptylem σ 2 pro všecha i (vizte poznámku 1 na str. 11). Cílem aproximace je najít funkci, jejíž graf „co nejlépe vystihuje rozložení těcho bodů (nemusí jimi procházet). Poté lze pomocí nalezeného funkčního předpisu počítat hodnoty pro libovolné x, samozřejmě s jistou chybou. Názor na to, co považovat za „co nejlepší vystižení rozložení bodů není jednoznačný. Kvůli poměrné jednoduchosti se nejčastěji používá metoda nejmenších čtverců (čtverec = druhá mocnina): Předpokládejme, že nejlépe vystihuje rozložení daných bodů [xi , yi ], i = 1, . . . , n, funkce y = f (x), jejíž funkční předpis obsahuje neznámé parametry. Pro každé i vypočteme čtverec (yi − f (xi ))2 . Sečteme všechny tyto čtverce, tj. vypočteme Φ=
n
(yi − f (xi ))2
i=1
a hledáme, pro které hodnoty parametrů je tato suma minimální, tj. kde je parciální derivace sumy Φ podle každého parametru 0. Pokusme se body [xi , yi ], i = 1, . . . , n, proložit přímku f (x) = β0 + β1 x, nazývanou též regresní.
4
Budeme tedy minimalizovat výraz n
Φ=
(yi − β1 xi − β0 )2 .
i=1
Vypočteme obě parciální derivace a upravíme je: ∂Φ = 2 (yi − β1 xi − β0 ) (−xi ) = 0 ∂β1 i=1 n
∂Φ = 2 (yi − β1 xi − β0 ) (−1) = 0 ∂β0 i=1 n
n
−xi yi + β1 x2i + β0 xi = 0
i=1 n
(−yi + β1 xi + β0 ) = 0
i=1 n
xi yi = β0
i=1
n
xi + β1
n
i=1 n
i=1
yi = nβ0 + β1
i=1
x2i
n
xi
i=1
n
n 2 x − x i i i=1 i=1 xi yi βˆ0 =
i=1
2 n n 2 n i=1 xi − ( i=1 xi )
n
n
n n x y − x i i i i=1 i=1 i=1 yi βˆ1 =
n 2
n 2 n i=1 xi − ( i=1 xi )
n
i=1 yi
n
(vypočtené koeficienty – odhady parametrů – je zvykem značit se stříškou) Získanou soustavu lze zapsat pomocí matic: ⎞ ⎛ ⎞ ⎛ n n yi ⎟ ⎜ n xi ⎟ ⎜ β0 ⎟ ⎜ i=1 ⎟ ⎜ i=1 × ⎟ = ⎜ ⎟ ⎜ n n n ⎠ β1 ⎠ ⎝ ⎝ xi yi xi x2i i=1
i=1
5
i=1
Tato soustava se dá sestavit pomocí těchto matic: ⎛ ⎞ ⎛ ⎞ y1 1 x1 ⎜ y2 ⎟ ⎜1 x2 ⎟ β0 ⎜ ⎟ ⎜ ⎟ X = ⎜ .. .. ⎟ , Y = ⎜ .. ⎟ , β = β1 ⎝.⎠ ⎝. . ⎠ 1 xn yn Soustava pak má tvar a
X Y = X Xβ β = (X X)−1 X Y.
Pokud bychom neprokládali přímku, že příslušné matice budou mít tvar ⎛ 1 x1 x21 . . . ⎜1 x2 x2 . . . 2 ⎜ X = ⎜ .. .. .. .. ⎝. . . .
ale polynom stupně n, projeví se to tak, ⎞ xn1 xn1 ⎟ ⎟ .. ⎟ , . ⎠
1 xn x2n . . . xnn
⎛
⎞ β0 ⎜β ⎟ ⎜ 1⎟ β = ⎜ . ⎟, ⎝ .. ⎠ βn
a opět β = (X X)−1 X Y. Aproximující křivka f (x) nemusí procházet danými body [xi , yi ], i = 1, . . . , n, je tedy rozumné spočítat odchylky ei = yi − f (xi ), nazývané rezidua. Dále se pak zkoumá součet čtverců reziduí Se = ni=1 e2i neboli reziduální součet čtverců. Platí, že rozptyl S 2 = Se /(n − k), kde k je počet parametrů (v našem případě 2), S je směrodatná odchylka. Kvalitu s jakou jsou data popsána regresní přímkou, udává index determinace n (ˆ yi − y¯)2 I = i=1 , n ¯)2 i=1 (yi − y kde yˆi je hodnota vypočtená z nalezené funkce, y¯ značí aritmetický průměr. Někdy se při výpočtu s výhodou používá také vztah y¯ = β0 + β1 x¯, který získáme sečtením rovnic yi = β0 + β1 xi pro všechna i a vydělení součtu n. Příklad 1. Mravenec průzkumník se probouzí při teplotě okolo 5 ◦ C, při teplotě 10 ◦ C už může dosáhnout rychlosti 18 m/hod., při teplotě 15 ◦ C vyvine rychlost 54 m/hod., při teplotě 20 ◦ C běží rychlostí 126 m/hod., při teplotě 25 ◦ C uhání rychlostí 210 m/hod., při teplotě 28 ◦ C jeho rychlost klesá na 190 m/hod. [Bernard Werber: Mravenci, KK, 2005] 1. Najděte regresní přímku β0 + β1 x pro závislost rychlosti y mravence průzkumníka na teplotě okolí x. Určete index determinace. 6
Obr. 3: Dané hodnoty rychlosti v závislosti na teplotě Řešení: Sestavíme tabulku i 1 2 3 4 5
xi 10 15 20 25 28 98
yi x2i xi yi 18 100 180 54 225 810 126 400 2520 210 625 5250 190 784 5320 598 2134 14080 βˆ0 =
=
n
i=1
yi2 yˆi ei e2i 324 13.37 4.630 21.441 2916 68.70 −14.698 216.030 15876 124.03 1.974 3.896 44100 179.35 30.645 939.140 36100 212.55 −22.552 508.570 99316 598 0 1689.1
n
x2i −
i=1 n ni=1 x2i −
yi
n i=1 xi i=1
2 ( ni=1 xi )
n
xi yi
=
−103708 598 · 2134 − 98 · 14080 = −97.287 = 2 5 · 2134 − 98 1066
n
n
n x y − x i i i i=1 i=1 yi βˆ1 = =
n 2
n 2 n i=1 xi − ( i=1 xi ) 11796 5 · 14080 − 98 · 598 = 11.066 = = 5 · 2134 − 982 1066 n
i=1
Směrodatná odchylka (vizte opět poznámku 1 na str. 11) je 1689.1 = 23.728. S= 3 7
Rovnice přímky je β0 + β1 x = −97.287 + 11.066x. Pro zjištění indexu determinace výpočteme průměry x ¯ = 19.6, y¯ = 119.6 a sestavíme další tabulku i 1 2 3 4 5
xi 10 15 20 25 28 98
yi yi − y¯ (yi − y¯)2 yˆi − y¯ (ˆ yi − y¯)2 18 −101.6 10322.56 −106.23 11284.81 54 −65.6 4303.36 −50.90 2590.81 126 6.4 40.96 4.43 19.26 210 90.4 8172.16 59.75 3570.06 190 70.4 4956.16 92.95 8639.70 598 0 27795.2 0 26105.00
Odtud I=
n (ˆ yi
i=1 n i=1 (yi
− y¯)2 = − y¯)2
26105.0 = 0.969 27795.2
Obr. 4: Regresní přímka 2. Najděte interpolační polynom rychlosti pohybu mravence v závislosti na teplotě. Řešení: Použijeme Lagrangeovu interpolaci pro daná data: teplota := [10, 15, 20, 25, 28], rychlost := [18, 54, 126, 210, 190]. 8
Vypočteme Lagrangeovy bázové polynomy a sestavíme polynom L(x): (x − 15)(x − 20)(x − 25)(x − 28) = (10 − 15)(10 − 20)(10 − 25)(10 − 28) = 0.0000740741x4 − 0.00651852x3 + 0.211481x2 − 2.99259x + 15.5556 (x − 10)(x − 20)(x − 25)(x − 28) = 1 (x) = (15 − 10)(15 − 20)(15 − 25)(15 − 28) = −0.000307692x4 + 0.0255385x3 − 0.766154x2 + 9.72308x − 43.0769 0 (x) =
(x − 10)(x − 15)(x − 25)(x − 28) = (20 − 10)(20 − 15)(20 − 25)(20 − 28) = 0.0005x4 − 0.039x3 + 1.0875x2 − 12.725x + 52.5
2 (x) =
(x − 10)(x − 15)(x − 20)(x − 28) = (25 − 10)(25 − 25)(25 − 20)(25 − 28) = −0.000444444x4 + 0.0324444x3 − 0.848889x2 + 9.42222x − 37.3333 (x − 10)(x − 15)(x − 20)(x − 25) = 4 (x) = (28 − 10)(28 − 15)(28 − 20)(28 − 25) = 0.000178063x4 − 0.0124644x3 + 0.316061x2 − 3.42771x + 13.3547
3 (x) =
L(x) = 180 (x) + 541(x) + 1262 (x) + 2103 (x) + 1904 (x) = = −0.0117835x4 + 0.792843x3 − 18.7557x2 + 195.232x − 733.761
200
150
100
50
10
15
20
25
t
Obr. 5: Interpolační polynom 9
30
Na zvěr ještě nakreslíme interpolační splajn rychlosti pohybu mravence v závislosti na teplotě získaný Matlabem. Kód úlohy byl tento: mravenci=[5,0;10,18;15,54;20,126;25,210;28,190]; X=mravenci(:,1); Y=mravenci(:,2); rad=length(X); xx=X(1):0.2:X(rad); yy=spline(X,Y,xx); plot(X,Y,’o’,xx,yy) Výsledek: 250
200
150
100
50
0
5
10
15
20
Obr. 6: Interpolační splajn
10
25
30
Poznámka 1. Význam směrodatné odchylky je znázorněn na obr. 7
Obr. 7: Křivka normálního rozdělení četnosti Procenta udávají plochu nad příslušným intervalem. Pravděpodobnost toho, že realizace náhodné veličiny s rozdělením N(μ, σ 2 ) bude v intervalu μ − σ, μ + σ je 68.27 %, v intervalu μ − 2σ, μ + 2σ je 95.45 % a v intervalu μ − 3σ, μ + 3σ je 99.73 %. Např. vyrábíme-li metrové tyče s přesností ±1 mm a zákazníkovi vyhovuje tyč délky 99 mm až 101 mm, pak vyrobíme pravděpodobně 31.73 % zmetků. Chceme-li méně zmetků, pak musíme zpřesnit výrobní proces, při přesnosti ±0.5 mm bude jen 4.55 % zmetků, při přesnosti ±0.33 mm bude jen 0.27 % zmetků. Také lze přesvědčit zákazníka, že na nějakém milimetru nesejde.
11
