5.1. Interpolaˇ cn´ı polynom

(1)

5. Interpolace a aproximace funkc´ı

Pr˚uvodce studiem Casto je potˇrebaˇ

”sloˇzitou” funkci f nahradit funkc´ı

”jednoduˇsˇs´ı”. V této kapitole budeme pˇredpokládat, ˇze u funkce f zn´ame jej´ı funkˇcn´ı hodnoty f_i = f (x_i) v uzlech x_i pro i = 0, . . . , n. Budeme rozliˇsovat dvˇe úlohy.

Interpolaˇcn´ı ´uloha: Hled´ame funkci ϕ, pro niˇz je

ϕ(x_i) = f_i, i = 0, . . . , n. (5.0.1) Aproximace metodou nejmenˇs´ıch ˇctverc˚u: Hled´ame funkci ϕ, pro niˇz je

ϕ(x_i)≈ fi, i = 0, . . . , n, (5.0.2) kde pˇribliˇzn´a rovnost

”≈” je urˇcena tak, aby souˇcet druchých mocnin odchylek mezi pˇredepsan´ymi hodnotami f_i a pˇredpokládan´ymi hodnotami ϕ(x_i) byl mi- nimáln´ı.

Jestliˇze tyto úlohy znázorn´ıme graficky, bude ˇreˇsen´ı interpolaˇcn´ı úlohy proch´azet pˇres body (x_i, f_i), i = 0, . . . , n, kdeˇzto ˇreˇsen´ı aproximaˇcn´ı úlohy bude (obecnˇe) procházet jejich bl´ızkým okol´ım.

Formulace obou úloh je zat´ım pˇr´ıliˇs obecná, protoˇze jsme neˇrekli jakého typu má b´yt funkce ϕ. Uk´aˇzeme tˇri volby: polynom, splajn (spline-funkce) a lineárn´ı kombinace obecných funkc´ı. Polynom je jednoduchý z hlediska matematických operac´ı (snadno se derivuje, integruje atp.), jeho graf vˇsak ˇcasto osciluje. Lepˇs´ı tvary grafu maj´ı splajny. Kombinace obecných funkc´ı se pouˇz´ıvá zpravidla v situac´ıch, kdy je známo, jakou závislost daná data popisuj´ı (pro periodickou závislost je dobré pouˇz´ıt funkce goniometrické, pro strmˇe rostouc´ı data se hod´ı funkce exponenciáln´ı atp.).

(2)

5.1. Interpolaˇ cn´ı polynom

C´ıle

Uk´aˇzeme metody pro sestaven´ı interpolaˇcn´ıho polynomu a odvod´ıme vzorec pro interpolaˇcn´ı chybu.

Pˇredpokl´adan´e znalosti

Polynomy. ˇReˇsen´ı soustav line´arn´ıch rovnic. Vˇeta o stˇredn´ı hodnotˇe dife- renci´aln´ıho poˇctu.

V´yklad

Funkci ϕ v ´uloze (5.0.1) budeme hledat jako interpolaˇcn´ı polynom stupnˇe nejv´yˇse n, tj. poloˇz´ıme ϕ = pn, kde

p_n(x) = a₀+ a₁x + a₂x²+ . . . + a_nxⁿ. (5.1.1) Zaˇcneme pˇr´ıkladem.

Pˇr´ıklad 5.1.1. Jsou d´any uzly x₀ = −2, x1 = −1, x2 = 1, x₃ = 2 a funkˇcn´ı hodnoty f₀ = 10, f₁ = 4, f₂ = 6, f₃ = 3. Urˇcete interpolaˇcn´ı polynom p₃.

Reˇˇ sen´ı: Hledaný polynom má obecný tvar

p₃(x) = a₀+ a₁x + a₂x² + a₃x³.

Koeﬁcienty a₀, a₁, a₂, a₃ urˇc´ıme tak, aby platilo (5.0.1). Kaˇzd´a interpolaˇcn´ı rovnost urˇcuje jednu rovnici:

p₃(−2) = 10 ⇒ a₀ − 2a₁ + 4a₂ − 8a₃ = 10, p₃(−1) = 4 ⇒ a₀ − a₁ + a₂ − a₃ = 4, p₃(1) = 6 ⇒ a0 + a₁ + a₂ + a₃ = 6, p₃(2) = 3 ⇒ a0 + 2a₁ + 4a₂ + 8a₃ = 3.

(3)

Dostali jsme soustavu line´arn´ıch rovnic

⎛

⎜⎜

⎜⎝

1 −2 4 −8 1 −1 1 −1

1 1 1 1

1 2 4 8

⎞

⎟⎟

⎟⎠

⎛

⎜⎜

⎜⎝ a₀ a₁ a₂ a₃

⎞

⎟⎟

⎟⎠

=

⎛

⎜⎜

⎜⎝ 10

4 6 3

⎞

⎟⎟

⎟⎠ ,

jej´ımˇz ˇreˇsen´ım (na tˇri desetinn´a m´ısta) jsou koeﬁcienty a₀ = 4.500, a₁ = 1.917, a₂ = 0.500 a a₃ =−0.917. Interpolaˇcn´ı polynom m´a tvar

p₃(x) = 4.500 + 1.917x + 0.500x²− 0.917x³. Jeho graf je na obr´azku 5.1.1.

−2 −1 0 1 2

2 4 6 8 10 12

Obr´azek 5.1.1: Graf interpolaˇcn´ıho polynomu p₃.

Rozborem postupu z pˇr´ıkaldu dok´aˇzeme n´asleduj´ıc´ı vˇetu.

Vˇeta 5.1.1.

Necht’ jsou d´any vz´ajemnˇe r˚uzn´e uzly x_i a funkˇcn´ı hodnoty f_i, i = 0, . . . , n.

Existuje pr´avˇe jeden interpolaˇcn´ı polynom stupnˇe nejv´yˇse n.

D˚ukaz: Dosazen´ım obecn´eho tvaru polynomu (5.1.1) do interpolaˇcn´ıch rovnost´ı

(4)

(5.0.1) dostaneme soustavu line´arn´ıch rovnic

p_n(x_i) = a₀+ a₁x_i + a₂x²_i + . . . + a_nxⁿ_i = f_i, i = 0, . . . , n, kterou lze zapsat maticovˇe jako

⎛

⎜⎜

⎝

1 x₀ x²₀ . . . xⁿ₀ 1 x₁ x²₁ . . . xⁿ₁ 1 x₂ x²₂ . . . xⁿ₂ ... ... ... . .. ...

1 x_n x²_n . . . xⁿ_n

⎞

⎟⎟

⎠

⎛

⎜⎜

⎝ a₀ a₁ a₂ ... a_n

⎞

⎟⎟

⎠

=

⎛

⎜⎜

⎝ f₀ f₁ f₂ ... f_n

⎞

⎟⎟

⎠ .

Matice této soustavy má nenulový determinant (Vandermod˚uv determinant).

Odtud plyne existence jediného ˇreˇsen´ı soustavy lineárn´ıch rovnic a také inter-

polaˇcn´ıho polynomu. 2

5.1.1. Lagrange˚uv tvar interpolaˇcn´ıho polynomu

Uk´aˇzeme postup, pˇri nˇemˇz se obejdeme bez ˇreˇsen´ı soustavy line´arn´ıch rovnic.

Interpolaˇcn´ı polynom budeme hledat ve tvaru

p_n(x) = f₀ϕ₀(x) + f₁ϕ₁(x) + . . . + f_nϕ_n(x). (5.1.2) Rovnosti p_n(x_i) = f_i, i = 0, 1, . . . , n budou splnˇeny, jestliˇze bude platit

ϕi(xj) =

⎧⎪

⎨

⎪⎩

1 pro i = j, 0 pro i= j.

Z vˇety 5.1.1. v´ıme, ˇze interpolaˇcn´ı polonom je stupnˇe nejv´yˇse n, takˇze také vˇsechny funkce ϕ_i mus´ı být polynomy stupnˇe nejv´yˇse n. Uveden´ym poˇzadavk˚um vyhovuje následuj´ıc´ı definice:

ϕ_i(x) = (x− x0) . . . (x− xi−1)(x− xi+1) . . . (x− xn)

(x_i− x₀) . . . (x_i− xi−1)(x_i − xi+1) . . . (x_i− xn) (5.1.3)

(5)

pro i = 0, 1, . . . , n. ˇCitatel je totiˇz polynom, který nabývá nulových hodnot ve vˇsech uzlech kromˇe x_i. V uzlu x_i pak nabývá nenulové hodnoty, která je obsaˇzena ve jmenovateli zlomku, takˇze plat´ı ϕ_i(x_i) = 1.

Polynom˚um ϕi, i = 0, 1, . . . , n se ˇr´ıká Lagrangeova báze interpolaˇcn´ı úlohy a vzorec (5.1.2) se nazýv´a Lagrange˚uv tvar inteprolaˇcn´ıho polynomu.

Pˇr´ıklad 5.1.2. Mˇejme d´any uzly x₀ =−2, x1 =−1, x2 = 1, x₃ = 2 a funkˇcn´ı hodnoty f₀ = 10, f₁ = 4, f₂ = 6, f₃ = 3. Napiˇste Lagrange˚uv tvar interpolaˇcn´ıho polynomu.

Reˇˇ sen´ı: Nejdˇr´ıve sestav´ıme Lagrangeovu b´azi. Podle (5.1.3) je

ϕ₀(x) = (x + 1)(x− 1)(x − 2)

(−2 + 1)(−2 − 1)(−2 − 2) = − 1

12(x + 1)(x− 1)(x − 2), ϕ₁(x) = (x + 2)(x− 1)(x − 2)

(−1 + 2)(−1 − 1)(−1 − 2) = 1

6(x + 2)(x− 1)(x − 2), ϕ₂(x) = (x + 2)(x + 1)(x− 2)

(1 + 2)(1 + 1)(1− 2) = −1

6(x + 2)(x + 1)(x− 2), ϕ₃(x) = (x + 2)(x + 1)(x− 1)

(2 + 2)(2 + 1)(2− 1) = 1

12(x + 2)(x + 1)(x− 1).

Dosazen´ım do (5.1.2) dostaneme v´ysledek

p₃(x) = −5

6(x + 1)(x− 1)(x − 2) +2

3(x + 2)(x− 1)(x − 2) −

−(x + 2)(x + 1)(x − 2) +1

4(x + 2)(x + 1)(x− 1).

Pozn´amka

Interpolaˇcn´ı polynom je podle vˇety 5.1.1. urˇcen jednoznaˇcnˇe. ´Upravou Lagran- geova tvaru proto mus´ıme nutnˇe doj´ıt k polynomu, kter´y jsem vypoˇc´ıtali v pˇr´ıkladu 5.1.1. (ovˇeˇrte).

(6)

5.1.2. Newton˚uv tvar interpolaˇcn´ıho polynomu Uvaˇzujme z´apis polynomu ve tvaru:

p_n(x) = a₀+a₁(x−x₀)+a₂(x−x₀)(x−x₁)+. . .+a_n(x−x₀) . . . (x−xn−1). (5.1.4) Jestliˇze dosad´ıme do interpolaˇcn´ıch rovnost´ı p_n(x_i) = f_i, i = 0, 1, . . . , n, dostaneme soustavu line´arn´ıch rovnic s doln´ı troj´uheln´ıkovou matic´ı:

⎛

⎜⎜

⎜⎝

1 0 0 . . . 0

1 (x₁− x₀) 0 . . . 0

1 (x₂− x0) (x₂− x0)(x₂− x1) . . . 0

... ... ... . .. ...

⎞

⎟⎟

⎟⎠

⎛

⎜⎜

⎜⎝ a₀ a₁ a₂ ...

⎞

⎟⎟

⎟⎠

=

⎛

⎜⎜

⎜⎝ f₀ f₁ f₂ ...

⎞

⎟⎟

⎟⎠

, (5.1.5)

Odud m˚uˇzeme postupnˇe vyj´adˇrit koeﬁcienty a_k: a₀ = f₀, a₁ = f₁ − a₀

x₁− x0 = f₁− f₀ x₁− x0,

a₂ = f₂− a₁(x₂− x₀)− a₀ (x₂− x0)(x₂− x1) =

f₂− f₁

x₂− x1 − f₁− f₀ x₁− x0

x₂ − x0 , atd..

Výrazy na pravých stranách jsou pomˇerné diference, jejichˇz oznaˇcen´ı zavád´ıme v následuj´ıc´ı definici.

Definice 5.1.1.

Necht’ jsou d´any vz´ajemnˇe r˚uzn´e uzly x_i a funkˇcn´ı hodnoty f_i, i = 0, . . . , n.

Pomˇern´e diference k-t´eho ˇr´adu f [x_i_+k, . . . , x_i], i = 0, 1, . . . , n− k deﬁnujeme rekurentnˇe:

• pro k = 0 : f [x_i] = f_i;

• pro k = 1 : f [x_i₊₁, x_i] = f_i₊₁− fi

x_i₊₁− xi

;

• pro k ≤ n : f[xi+k, . . . , , x_i] = f [x_i_+k, . . . , x_i₊₁]− f[xi+k−1, . . . , x_i] x_i_+k− xi

.

(7)

Porovn´an´ım pomˇern´ych diferenc´ı s koeﬁcienty a_k vid´ıme, ˇze a_k = f [x_k, . . . , x₀], k = 0, 1, . . . , n.

Dosazen´ım do (5.1.4) dostaneme Newton˚uv tvar interpolaˇcn´ıho polynomu:

p_n(x) = f₀+ f [x₁, x₀](x− x₀) + . . . + f [x_n, . . . , x₀](x− x₀) . . . (x− xn−1). (5.1.6) Pˇri jeho sestavován´ı potˇrebujeme vypoˇc´ıtat pomˇerné diference. Vˇse ukáˇzeme v ná- sleduj´ıc´ım pˇr´ıkladu.

Pˇr´ıklad 5.1.3. Mˇejme d´any uzly x₀ =−2, x₁ =−1, x₂ = 1, x₃ = 2 a funkˇcn´ı hodnoty f₀ = 10, f₁ = 4, f₂ = 6, f₃ = 3. Napiˇste Newton˚uv tvar interpolaˇcn´ıho polynomu.

Reˇˇ sen´ı: Potˇrebujeme vypoˇc´ıtat pomˇern´e diference:

f [x₁, x₀], f [x₂, x₁, x₀], f [x₃, x₂, x₁, x₀].

Podle deﬁnice je

f [x₁, x₀] = 4− 10

−1 + 2 = −6, f [x₂, x₁, x₀] = f [x₂, x₁]− f[x1, x₀]

x₂− x₀ =

6−41+1+ 6 1 + 2 = 7

3,

f [x₃, x₂, x₁, x₀] = f [x₃, x₂, x₁]− f[x2, x₁, x₀]

x₃− x0 =

3−62−1−⁶⁻⁴₁₊₁ 2+1 − ⁷₃

2 + 2 = −11 12. Dosazen´ım do (5.1.6) dostaneme v´ysledek

p₃(x) = 10− 6(x + 2) + 7

3(x + 2)(x + 1)−11

12(x + 2)(x + 1)(x− 1).

Pˇrehlednˇe m˚uˇzeme výpoˇcet pomˇerných diferenc´ı provést v tabulce (tabulka 5.1.1), kde do prvn´ıch dvou sloupc˚u zap´ıˇseme zadané uzly a funkˇcn´ı hodnoty a v kaˇzdém dalˇs´ım sloupci pak vypoˇc´ıtáme vˇsechny (!) pomˇerné diference postupnˇe se zvy- ˇsuj´ıc´ıch ˇrád˚u. Pro napsán´ı interpolaˇcn´ıho polynomu potˇrebujeme z této tabulky hodnoty diferenc´ı z prvn´ıho ˇrádku.

(8)

Tabulka 5.1.1: V´ypoˇcet pomˇern´ych diferenc´ı.

i x_i f_i f [x_i₊₁, x_i] f [x_i₊₂, x_i₊₁, x_i] f [x₃, x₂, x₁, x₀]

0 −2 10 −6 ⁷₃ −¹¹₁₂

1 −1 4 1 −⁴₃

2 1 6 −3

3 2 3

5.1.3. Interpolaˇcn´ı chyba

Pˇredpokl´adejme, ˇze hodnoty f_i jsou funkˇcn´ımi hodnotami funkce f v uzlech x_i, tj. f_i = f (x_i). Bude n´as zaj´ımat interpolaˇcn´ı chyba

f (x)− pn(x).

V uzlech x_i je interpolaˇcn´ı chyba nulová, ale mimo uzly m˚uˇze být velká.

Vˇeta 5.1.2.

Necht’ uzly x_i, i = 0, 1, . . . , n, jsou vz´ajemnˇe r˚uzn´e a leˇz´ı na intervalua, b. Necht’

funkce f m´a na tomto intervalu n + 1 spojit´ych derivac´ı. Pak pro kaˇzd´e x∈ a, b

existuje ξ = ξ(x) v (a, b) tak, ˇze plat´ı

f (x)− pn(x) = f⁽ⁿ⁺¹⁾(ξ)

(n + 1)! π_n₊₁(x), (5.1.7) kde π_n₊₁(x) = (x− x0) . . . (x− xn).

D˚ukaz: Pro x = x_i je rovnost (5.1.7) splnˇena, protoˇze obˇe jej´ı strany jsou nulov´e.

Pro pevnˇe zvolen´e x= xi deﬁnujme funkci

g(t) = f (t)− pn(t)− π_n₊₁(t)

π_n₊₁(x)(f (x)− pn(x)) , (5.1.8) kde t je promˇenn´a a x je parametr. Funkce g m´a zˇrejmˇe n+2 koˇren˚u, kterými jsou body x₀, . . ., x_na x. Kaˇzd´a derivace funkce g m´a o jeden koˇren ménˇe, takˇze (n+1)- n´ı derivace má jediný koˇren v nˇejakém bodˇe ξ ∈ (a, b). Derivujeme-li (n + 1)-krát

(9)

v´yraz (5.1.8) (podle t) a pouˇzijeme pˇritom p⁽ⁿ⁺¹⁾n (t) = 0 a π⁽ⁿ⁺¹⁾_n₊₁ (t) = (n + 1)!, dostaneme

0 = g⁽ⁿ⁺¹⁾(ξ) = f⁽ⁿ⁺¹⁾(ξ)− (n + 1)!

π_n₊₁(x)(f (x)− pn(x)) .

Jestliˇze odtud vyjádˇr´ıme interpolaˇcn´ı chybu, vznikne rovnost (5.1.7). 2 Na pr˚ubˇeh interpolaˇcn´ı chyby v intervalu a, b má podstatný vliv tvar poly- nomu π_n₊₁, jak ukazuje následuj´ıc´ı pˇr´ıklad.

Pˇr´ıklad 5.1.4. (Rungeho pˇr´ıklad) Nakresl´ıme graf funkce f (x) = 1

1 + x²

a graf interpolaˇcn´ıho polynomu odpov´ıdaj´ıc´ıho uzl˚um xi =−5+i, i = 0, 1, . . . , 10.

V´ysledek porovn´ame s grafem polynomu

π₁₁(x) = (x + 5)(x + 4) . . . (x− 5).

Reˇˇ sen´ı: Obr´azek 5.1.2.a ukazuje graf polynomu π₁₁. Z jeho pr˚ubˇehu lze usou- dit, ˇze nejvˇetˇs´ı interpolaˇcn´ı chyby budou pobl´ıˇz krajn´ıch uzl˚u x₀ =−5 a x₁₀ = 5.

Na obr´azku 5.1.2.b vid´ıme, ˇze graf interpolaˇcn´ıho polynomu osciluje kolem grafu funkce f a ˇze oscilace jsou nejvˇetˇs´ı pr´avˇe na kraj´ıch intervalu −5, 5. Pozna- menejme jeˇstˇe, ˇze pˇri zvˇetˇsen´ı poˇctu interpolaˇcn´ıch uzl˚u nedojde ke zmenˇsen´ı interpolaˇcn´ı chyby, ale naopak k jej´ımu zvˇetˇsen´ı.

Kontroln´ı ot´azky

Otázka 1. Jaké znáte metody pro sestaven´ı interpolaˇcn´ıho polynomu?

Ot´azka 2. Jak´eho stupnˇe je interpolaˇcn´ı polynom?

Ot´azka 3. Jak se chov´a interpolaˇcn´ı chyba?

Ulohy k samostatn´´ emu ˇreˇsen´ı

1. Pro uzly x₀ =−1, x1 = 0, x₂ = 2, x₃ = 3, x₄ = 5 a funkˇcn´ı hodnoty f₀ =−2,

(10)

−5 0 5

−5 0 5x 10⁵

−5 0 5

−0.5 0 0.5 1 1.5 2

a b

Obr´azek 5.1.2: a) Graf π₁₁; b) Grafy f (neosciluj´ıc´ı) a p₁₀ (osciluj´ıc´ı).

f₁ = 1, f₂ = 0, f₃ = 2, f₄ =−1 vypoˇctˇete interpolaˇcn´ı polynom ve tvaru (5.1.1).

2. Pro pˇredchoz´ı data vypoˇctˇete Lagrange˚uv a Newton˚uv tvar interpolaˇcn´ıho polynomu.

Výsledky úloh k samostatnému ˇreˇsen´ı 1. p₄(x) =−₂₀³ x⁴+¹¹₁₀x³ −¹⁰⁹₆₀x²− ₁₅¹ x + 1.

2. Lagrange˚uv tvar: p₄(x) =−₃₆¹x(x− 2)(x − 3)(x − 5) −₃₀¹(x + 1)(x− 2)(x − 3) (x− 5) − ₁₂¹(x + 1)x(x− 2)(x − 5) − ₁₈₀¹ (x + 1)x(x− 2)(x − 3);

Newton˚uv tvar: p₄(x) =−2 + 3(x + 1) −³⁵₃₀(x + 1)x +¹₂(x + 1)x(x− 2) −₂₀³ (x + 1) x(x− 2)(x − 3).

(11)

5.2. Interpolaˇ cn´ı splajny

C´ıle

Vidˇeli jsme, ˇze graf interpolaˇcn´ıho polynomu m˚uˇze nepˇr´ıjemnˇe oscilovat. Tato situace nastává pˇri pˇredepsán´ı vˇetˇs´ıho poˇctu dat, protoˇze interpolaˇcn´ı polynom je pak vysokého stupnˇe. Zdá se proto rozumné pˇri ˇreˇsen´ı interpolaˇcn´ı úlohy pouˇz´ıt funkci, kter´a bude poˇcástech polynomem n´ızk´eho stupnˇe, jej´ıˇz jednotlivé ˇcásti budou na sebe navazovat dostateˇcnˇe hladce. Takovým funkc´ım se ˇr´ık´a splajn (z angl.

”spline”). Ukáˇzeme dva nejˇcastˇeji pouˇz´ıvané splajny: lineárn´ı a kubický.

Interpolaˇcn´ı polynom. Spojitost derivace. ˇReˇsen´ı soustav line´arn´ıch rovnic.

V´yklad

Abychom se vyhnuli komplikac´ım pˇri popisu, budeme pˇredpokl´adat, ˇze uzly interpolace tvoˇr´ı rostouc´ı posloupnost, tzn. x₀ < x₁ < . . . < x_n. Vzd´alenost dvou sousedn´ıch uzl˚u oznaˇc´ıme h_i, tj. h_i = x_i − xi−1, i = 1, . . . , n.

5.2.1. Line´arn´ı splajn Definice 5.2.1.

Lineárn´ım splajnem nazýv´ame funkci s₁, která je spojitá na intervalu x0, x_n a na kaˇzdém podintervaluxi−1, x_i, i = 1, . . . , n, je polynomem prvn´ıho stupnˇe.

Line´arn´ı interpolaˇcn´ı splajn je ˇreˇsen´ım ´ulohy (5.0.1), tzn. ˇze pro nˇej plat´ı s₁(x_i) = f_i, i = 0, . . . , n. M˚uˇzeme ho zapsat po ˇc´astech pro i = 1, . . . , n:

s₁(x) = f_i₋₁(1− t) + fit, t = (x− xi−1)/h_i, x∈ xi−1, x_i. (5.2.1) Grafem line´arn´ıho splajnu je lomen´a ˇc´ara.

(12)

Pˇr´ıklad 5.2.1. Mˇejme d´any uzly x₀ =−2, x1 =−1, x2 = 1, x₃ = 2 a funkˇcn´ı hodnoty f₀ = 10, f₁ = 4, f₂ = 6, f₃ = 3. Napiˇste line´arn´ı interpolaˇcn´ı splajn.

Reˇˇ sen´ı: Zap´ıˇseme jej pomoc´ı pˇredpisu (5.2.1):

s₁(x) =

⎧⎪

⎪⎪

⎪⎨

⎪⎪

⎩

10ϕ₁(t) + 4ϕ₂(t), t = x + 2 pro x∈ −2, −1, 4ϕ₁(t) + 6ϕ₂(t), t = (x + 1)/2 pro x∈ −1, 1, 6ϕ₁(t) + 3ϕ₂(t), t = x− 1 pro x∈ 1, 2, kde ϕ₁(t) = 1− t, ϕ2(t) = t. Graf je zn´azornˇen na obr´azku 5.2.1.

5.2.2. Kubick´y splajn Definice 5.2.2.

Kubickým splajnem nazýv´ame funkci s₃, která má na intervalu x₀, x_n dvˇe spojité derivace a na kaˇzdém podintervaluxi−1, x_i, i = 1, . . . , n, je polynomem tˇret´ıho stupnˇe.

Kubick´y interpolaˇcn´ı splajn, je ˇreˇsen´ı interpolaˇcn´ı úlohy (5.0.1). Jeho kon- strukce je sloˇzitˇejˇs´ı neˇz u lineárn´ıho splajnu. Vyjdeme opˇet z vyjádˇren´ı po ˇcástech pro i = 1, . . . , n:

s₃(x) = f_i₋₁(1− 3t²+ 2t³) + f_i(3t²− 2t³)

+m_i₋₁h_i(t− 2t²+ t³) + m_ih_i(−t²+ t³), (5.2.2) kde t = (x− xi−1)/h_i a x∈ xi−1, x_i. Tento pˇredpis je navrˇzen tak, aby parame- try f_i₋₁, f_i a m_i₋₁, m_i mˇely v´yznam funkˇcn´ıch hodnot a hodnot prvn´ı derivace v krajn´ıch bodech intervalu xi−1, x_i, tj. plat´ı

s₃(x_i₋₁) = f_i₋₁, s₃(x_i) = f_i, (5.2.3) s₃(x_i₋₁) = m_i₋₁, s₃(x_i) = m_i. (5.2.4)

(13)

O splnˇen´ı rovnost´ı (5.2.3) a (5.2.4) se m˚uˇzeme pˇresvˇedˇcit dosazen´ım x_i₋₁ a x_i do (5.2.2) a do prvn´ı derivace s₃, kterou vyjádˇr´ıme z (5.2.2) podle pravidla o deri- vován´ı sloˇzené funkce:

s₃(x) = fi−1(−6t + 6t²)/hi+ fi(6t− 6t²)/hi

+m_i₋₁(1− 4t + 3t²) + m_i(−2t + 3t²). (5.2.5) Pˇredpis (5.2.2) zaruˇcuje spojitost prvn´ı derivace s₃ na celém intervalu x0, x_n pro libovoln´e hodnoty m_i. Spojitost druhé derivace vynut´ıme speciáln´ı volbou m_i. Budeme poˇzadovat

x→xlimi−s₃(x) = lim

x→xi+s₃(x) (5.2.6) ve vnitˇrn´ıch uzlech x_i, i = 1, . . . , n− 1. Potˇrebný výraz pro druhou derivaci vypoˇcteme z (5.2.5) opˇet podle pravidla o derivován´ı sloˇzené funkce:

s₃(x) = f_i₋₁(−6 + 12t)/h²i + f_i(6− 12t)/h²i

+m_i₋₁(−4 + 6t)/hi+ m_i(−2 + 6t)/hi. (5.2.7) Levou stranu v (5.2.6) vyj´adˇr´ıme z (5.2.7) pro t = 1:

x→xlimi−s₃(x) = 6fi−1/h²_i − 6fi/h²_i + 2mi−1/hi+ 4mi/hi. (5.2.8) Pravou stranu v (5.2.6) vyj´adˇr´ıme z (5.2.7) pro t = 0, kdyˇz souˇcasnˇe posuneme indexov´an´ı:

x→xlimi+s(x) = −6fi/h²_i₊₁+ 6fi+1/h²_i₊₁− 4mi/hi+1− 2mi+1/hi+1. (5.2.9) Dosad´ıme-li (5.2.8) a (5.2.9) do (5.2.6), dostaneme po jednoduch´e ´upravˇe

h_i₊₁m_i₋₁+ 2(h_i₊₁+ h_i)m_i+ h_im_i₊₁ = 3 −h_i₊₁

h_i f_i₋₁+ h_i₊₁

h_i − h_i h_i₊₁

f_i+ h_i h_i₊₁f_i₊₁

!

, i = 1, . . . , n− 1. (5.2.10)

(14)

Tyto rovnosti tvoˇr´ı soustavu n−1 rovnic pro n+ 1 nezn´am´ych mi, i = 0, 1, . . . .n.

Abychom dostali jedin´e ˇreˇsen´ı, urˇc´ıme m₀ a m_nnapˇr´ıklad jako pˇribliˇzn´e derivace:

m₀ = f₁− f₀

h₁ , m_n= f_n− fn−1

h_n . (5.2.11)

Pˇr´ıklad 5.2.2. Mˇejme d´any uzly x₀ =−2, x1 =−1, x2 = 1, x₃ = 2 a funkˇcn´ı hodnoty f₀ = 10, f₁ = 4, f₂ = 6, f₃ = 3. Napiˇste kubick´y interpolaˇcn´ı splajn.

Reˇˇ sen´ı: Nejdˇr´ıve vypoˇc´ıt´ame parametry m_i, i = 0, 1, 2, 3. Podle (5.2.11) je m₀ = 4− 10

−1 + 2 =−6, m3 = 3− 6

2− 1 =−3.

Soustava (5.2.10) m´a dvˇe rovnice:

2(h₂+ h₁)m₁+ h₁m₂ = 3 −h₂ h₁f₀+

h₂ h₁ − h₁

h₂

f₁+h₁ h₂f₂

!

− h2m₀,

h₃m₁ + 2(h₃+ h₂)m₂ = 3 −h₃ h₂f₁+

h₃ h₂ − h₂

h₃

f₂+h₂ h₃f₃

!

− h2m₃, kter´e m˚uˇzeme ps´at jako

⎛

⎜⎝ 6 1 1 6

⎞

⎟⎠

⎛

⎜⎝ m₁ m₂

⎞

⎟⎠ =

⎛

⎜⎝ −21

−9

⎞

⎟⎠ .

Vyˇreˇsen´ım dostaneme m₁ = −²³⁴₇₀, m₂ = −⁶⁶₇₀. Výsledný splajn zap´ıˇseme podle (5.2.2) po ˇcástech:

s₃(x) =

⎧⎪

⎪⎪

⎪⎨

⎪⎪

⎩

10ϕ₁(t) + 4ϕ₂(t)− 6ϕ₃(t)− ¹¹⁷₃₅ϕ₄(t), t = x + 2 pro x∈ −2, −1, 4ϕ₁(t) + 6ϕ₂(t)− ²³⁴₃₅ϕ₃(t)− ⁶⁶₃₅ϕ₄(t), t = (x + 1)/2 pro x∈ −1, 1, 6ϕ₁(t) + 3ϕ₂(t)− ³³₃₅ϕ₃(t)− 3ϕ₄(t), t = x− 1 pro x ∈ 1, 2, kde ϕ₁(t) = 1− 3t²+ 2t³, ϕ₂(t) = 3t²− 2t³, ϕ₃(t) = t− 2t²+ t³, ϕ₄(t) =−t²+ t³. Graf je zn´azornˇen na obr´azku 5.2.1.

Pˇr´ıklad 5.2.3. (Rungeho pˇr´ıklad, pokraˇcov´an´ı) Nakresl´ıme graf inter- polaˇcn´ıho kubick´eho splajnu pro funkci f a uzly x_i z pˇr´ıkladu 5.1.4. a porovn´ame ho s grafem interpolaˇcn´ıho polynomu.

(15)

−2 −1 0 1 2 2

4 6 8 10

s₁ s3

Obr´azek 5.2.1: Graf line´arn´ıho (s₁) a kubick´eho interpolaˇcn´ıho (s₃) splajnu.

Reˇˇ sen´ı: Na obr´azku 5.2.2 vid´ıme, ˇze splajn s₃ neosciluje a je proto mno- hem lepˇs´ı aproximac´ı interpolovan´e funkce f neˇz interpolaˇcn´ı polynom p₁₀, viz obr´azek 5.1.2.b.

−5 0 5

0 0.2 0.4 0.6 0.8 1

−5 0 5

0 0.2 0.4 0.6 0.8 1

a b

Obr´azek 5.2.2: a) Funkce f ; b) Kubick´y interpolaˇcn´ı splajn s₃.

Otázka 1. Co je to splajn? Jak se definuje a poˇc´ıtá splajn lineárn´ı a kubický?

Ot´azka 2. Jak se chovaj´ı pˇri interpolaci splajny v porovn´an´ı s polynomy?

(16)

1. Pro uzly x₀ =−1, x1 = 0, x₂ = 2, x₃ = 3, x₄ = 5 a funkˇcn´ı hodnoty f₀ =−2, f₁ = 1, f₂ = 0, f₃ = 2, f₄ =−1 sestavte line´arn´ı interpolaˇcn´ı splajn.

2. Pro stejn´a data sestavte kubick´y interpolaˇcn´ı splajn.

Výsledky úloh k samostatnému ˇreˇsen´ı 1.

s₁(x) =

⎧⎪

⎪⎪

⎪⎨

⎪⎪

⎩

−2ϕ1(t) + ϕ₂(t), t = x + 1 pro x∈ −1, 0, ϕ₁(t), t = x/2 pro x∈ 0, 2, 2ϕ₂(t), t = x− 2 pro x∈ 2, 3, 2ϕ₁(t)− 1ϕ2(t), t = (x− 3)/2 pro x ∈ 3, 5, kde ϕ₁(t) = 1− t, ϕ₂(t) = t.

2. Krajn´ı parametry jsou m₀ = 3, m₄ =−³₂, ostatn´ı dostaneme ze soustavy

⎛

⎜⎜

⎝

6 1 0 1 6 2 0 2 6

⎞

⎟⎟

⎠

⎛

⎜⎜

⎝ m₁ m₂ m₃

⎞

⎟⎟

⎠=

⎛

⎜⎜

⎝

212 212

9

⎞

⎟⎟

⎠,

takˇze m₁ = ⁹⁷₆₂, m₂ = ⁶⁹₆₂, m₃ = ³⁵₃₁ a koneˇcnˇe

s₃(x) =

⎧⎪

⎪⎪

⎪⎨

⎪⎪

⎩

−2ϕ₁(t) + ϕ₂(t) + 3ϕ₃(t) + ⁹⁷₆₂ϕ₄(t), t = x + 1 pro x∈ −1, 0, ϕ₁(t) + ⁹⁷₃₁ϕ₃(t) +⁶⁹₃₁ϕ₄(t), t = x/2 pro x∈ 0, 2, 2ϕ₂(t) + ⁶⁹₆₂ϕ₃(t) +³⁵₃₁3ϕ₄(t), t = x− 2 pro x ∈ 2, 3, 2ϕ₁(t)− ϕ2(t) +⁷⁰₃₁ϕ₃(t)− 3ϕ4(t), t = (x− 3)/2 pro x ∈ 3, 5, kde ϕ₁(t) = 1−3t²+ 2t³, ϕ₂(t) = 3t²−2t³, ϕ₃(t) = t−2t²+ t³, ϕ₄(t) =−t²+ t³.

(17)

5.3. Aproximace metodou nejmenˇ s´ıch ˇ ctverc˚ u

C´ıle

V mnoha situac´ıch, v nichˇz je potˇreba danou funkci f nahradit funkc´ı

”jed- noduˇsˇs´ı”, je nevhodné nebo v˚ubec nelze pouˇz´ıt interpolaci. Jsou-li napˇr´ıklad v uzlech zadány nepˇresné hodnoty, pˇrenáˇs´ı se tato nepˇresnost i na interpolant. Inter- polace je nepouˇzitelná, jestliˇze je poˇzadován jistý charakter aproximuj´ıc´ı funkce a pˇritom ˇzádná funkce tohoto charakteru nen´ı interpolantem. V tˇechto pˇr´ıpadech je rozumné pouˇz´ıt metodu nejmenˇs´ıch ˇctverc˚u.

Lineárn´ı závislost a nezávislost. Urˇcen´ı minima funkce pomoc´ı derivace. ˇReˇsen´ı soustav lineárn´ıch rovnic.

V´yklad

Zaˇcneme pˇr´ıkladem.

Pˇr´ıklad 5.3.1. Mˇejme d´any uzly x₀ =−2, x1 =−1, x2 = 1, x₃ = 2 a funkˇcn´ı hodnoty f₀ = 10, f₁ = 4, f₂ = 6, f₃ = 3. Najdˇete pˇr´ımku

ϕ(x) = c₁ + c₂x, (5.3.1)

kter´a je

”bl´ızko” pˇredepsan´ym hodnot´am.

Reˇˇ sen´ı: Nejdˇr´ıve se mus´ıme rozhodnout jak ch´apat slovo

”bl´ızko”. Uˇz jsme to vlastnˇe ˇrekli, kdyˇz jsme popisovali smysl pˇribliˇzn´ych rovnost´ı v aproximaˇcn´ı

´

uloze (5.0.2). Pˇr´ımku ϕ urˇc´ıme tak, aby minimalizovala souˇcet druh´ych mocnin odchylek₃

i=0(ϕ(x_i)−fi)². Jestliˇze sem dosad´ıme (5.3.1), dostaneme ´ulohu na mi- nimalizaci funkce dvou promˇenn´ych Ψ(c₁, c₂) = ₃

i=0(c₁+ c₂x_i− fi)². Minimum

(18)

c^∗₁, c^∗₂ vyhovuje rovnic´ım

∂Ψ

∂c₁(c^∗₁, c^∗₂) = 0, ∂Ψ

∂c₂(c^∗₁, c^∗₂) = 0.

Po vyjádˇren´ı parciáln´ıch derivac´ı dostáváme

2

3 i=0

(c^∗₁+ c^∗₂xi− fi) = 0, 2

3 i=0

(c^∗₁+ c^∗₂xi − fi)xi = 0,

coˇz je soustava line´arn´ıch rovnic

⎛

⎜⎜

⎝

3 i=0

1

3 i=0

x_i

3 i=0

x_i

3 i=0

x²_i

⎞

⎟⎟

⎠

⎛

⎝ c^∗₁ c^∗₂

⎞

⎠ =

⎛

⎜⎜

⎝

3 i=0

f_i

3 i=0

f_ix_i

⎞

⎟⎟

⎠, tj.

⎛

⎝ 4 0 0 10

⎞

⎠

⎛

⎝ c^∗₁ c^∗₂

⎞

⎠ =

⎛

⎝ 23

−12

⎞

⎠ .

Tato soustava m´a jedin´e ˇreˇsen´ı c^∗₁ = ²³₄, c^∗₂ =−⁶₅, takˇze hledan´a pˇr´ımka ϕ^∗ = ϕ je urˇcena pˇredpisem

ϕ^∗(x) = 23 4 − 6

5x. (5.3.2)

Jej´ı graf je zn´azornˇen na obr´azku 5.3.1. 2

−2 −1 0 1 2

2 4 6 8 10

Obr´azek 5.3.1: Aproximace metodou nejmenˇs´ıch ˇctverc˚u; pˇr´ımka (5.3.2) plnˇe;

funkce (5.3.8) ˇc´arkovanˇe.

(19)

Postup z pˇr´ıkladu nyn´ı zobecn´ıme. Budeme pˇredpokládat, ˇze je dán systém funkc´ı ϕ_j = ϕ_j(x), j = 1, . . . , m a budeme uvaˇzovat vˇsechny funkce ve tvaru

ϕ(x) = c₁ϕ₁(x) + . . . + c_mϕ_m(x) =

m j=1

c_jϕ_j(x), (5.3.3)

kde koeﬁcienty c₁, . . . , c_m jsou libovoln´a ˇc´ısla. Funkci ϕ^∗, pro niˇz plat´ı

n i=1

(ϕ^∗(x_i)− fi)² ≤

n i=1

(ϕ(x_i)− fi)² ∀ϕ (5.3.4)

naz´yv´ame aproximac´ı podle metody nejmenˇs´ıch ˇctverc˚u. Jej´ı koeﬁcienty c^∗₁, . . . , c^∗_m urˇc´ıme jako minimum funkce

Ψ(c₁, . . . , c_m) =

n i=1

(

m j=1

c_jϕ_j(x_i)− fi)², (5.3.5)

kter´e vyhovuje rovnic´ım

∂Ψ

∂c_k(c^∗₁, . . . , c^∗_m) = 0, k = 1, . . . , m. (5.3.6) Vyj´adˇr´ıme-li parci´aln´ı derivace

∂Ψ

∂c_k = 2

n i=1

(

m j=1

cjϕj(xi)− fi)ϕk(xi)

a dosad´ıme je do (5.3.6), dostaneme po jednoduché úpravˇe soustavu lineárn´ıch rovnic

m j=1

_n

i=1

ϕ_j(x_i)ϕ_k(x_i)

c^∗_j =

n i=1

f_iϕ_k(x_i), k = 1, . . . , m. (5.3.7)

Soustava (5.3.7) se naz´yv´a soustava norm´aln´ıch rovnic.

Vˇeta 5.3.1.

Necht’ jsou d´any vz´ajemnˇe r˚uzn´e uzly xi a funkˇcn´ı hodnoty fi, i = 0, . . . , n.

Necht’ je dán syst´em funkc´ı ϕ_j, j = 1, . . . , m, které jsou lineárnˇe nezávislé. Potom existuje jedin´a funkce ϕ^∗, která splˇnuje (5.3.4) a jej´ı koeficienty c^∗₁, . . . , c^∗_m jsou ˇreˇsen´ım soustavy normáln´ıch rovnic (5.3.7).

(20)

D˚ukaz: V bodˇe c^∗₁, . . . , c^∗_m, který vyhovuje rovnic´ım (5.3.6), nabývá funkce Ψ minima, jestliˇze matice druhých derivac´ı je symmetrická a pozitivnˇe definitn´ı (kladná). Druhé derivace jsou urˇceny vzorcem

∂²Ψ

∂c_k∂c_l = 2

n i=1

ϕl(xi)ϕk(xi),

odkud je symetrie zˇrejm´a na prvn´ı pohled (prohozen´ım index˚u k a l se nic nezmˇen´ı). Necht’ d₁, . . . , d_m jsou ˇc´ısla ne vˇsechny souˇcasnˇe nulov´a. Potom

m k=1

m l=1

d_kd_l ∂²Ψ

∂c_k∂c_l = 2

n i=1

"^m

l=1

d_lϕ_l(x_i)

#"^m

k=1

d_kϕ_k(x_i)

#

= 2

n i=1

˜

ϕ(x_i)² > 0, kde ˜ϕ(x) = _m

k=1d_kϕ_k(x), takˇze matice druh´ych derivac´ı je pozitivnˇe deﬁnitn´ı.

Odtud také plyne, ˇze matice soustavy normáln´ıch rovnic je regulárn´ı, coˇz zna- mená, ˇze existuje jej´ı jedin´e ˇreˇsen´ı c^∗₁, . . . , c^∗_m, které urˇcuje jedinou funkci ϕ^∗. 2

Pˇri aproximaci metodou nejmenˇs´ıch ˇctverc˚u se mus´ıme nejdˇr´ıve rozhodnout pro nˇejaký lineárnˇe nezávislý syst´em funkc´ı ϕ_j, j = 1, . . . , m. Poté staˇc´ı sestavit a vyˇreˇsit soustavu normáln´ıch rovnic (5.3.7).

Pˇr´ıklad 5.3.2. Napiˇste normáln´ı soustavu lineárn´ıch rovnic odpov´ıdaj´ıc´ı systému funkc´ı

ϕ₁(x) = e^−x, ϕ₂(x) = sin x.

Aproximujte data z pˇr´ıkladu 5.3.1.

Reˇˇ sen´ı: Obecnˇe m´a soustava norm´aln´ıch rovnic tvar

⎛

⎜⎜

⎜⎝

n i=0

e^−2xⁱ

n i=0

e^−xⁱsin x_i

n i=0

e^−xⁱsin x_i

n i=0

sin²x_i

⎞

⎟⎟

⎟⎠

⎛

⎝ c^∗₁ c^∗₂

⎞

⎠ =

⎛

⎜⎜

⎜⎝

n i=0

f_ie^−xⁱ

n i=0

f_isin x_i

⎞

⎟⎟

⎟⎠.

Po dosazen´ı dostaneme

⎛

⎜⎝ 62.1409 −8.5736

−8.5736 3.0698

⎞

⎟⎠

⎛

⎝ c^∗₁ c^∗₂

⎞

⎠ =

⎛

⎜⎝ 87.3770

−4.6821

⎞

⎟⎠

(21)

a odtud vypoˇc´ıt´ame c^∗₁ = 1.9452, c^∗₂ = 3.9076, tj.

ϕ^∗(x) = 1.9452e^−x+ 3.9076 sin x. (5.3.8) Graf je zn´azornˇen na obr´azku 5.3.1.

Ot´azka 1. Kdy je vhodn´e pouˇz´ıt metodu nejmenˇs´ıch ˇctverc˚u?

Otázka 2. Graficky znázornˇete smysl výrazu pro souˇcet druhých mocnin odchylek?

Otázka 3. Co je to normáln´ı soustava lineárn´ıch rovnic a jak vznikne?

Ot´azka 4. Co se stane, kdyˇz v (5.3.3) a (5.3.4) bude m = n + 1?

1. Napiˇste soustavu norm´aln´ıch line´arn´ıch rovnic pro syst´em funkc´ı ϕ₁(x) = 1, ϕ₂(x) = x, ϕ₃(x) = x².

2. Data x₀ = −1, x1 = 0, x₂ = 2, x₃ = 3, x₄ = 5 a f₀ = −2, f1 = 1, f₂ = 0, f₃ = 2, f₄ =−1 aproximujte metodou nejmenˇc´ıch ˇctverc˚u pomoc´ı syst´emu funkc´ı z pˇredchoz´ı ´ulohy.

Výsledky úloh k samostatnému ˇreˇsen´ı 1. Norm´aln´ı soustava má tvar:

⎛

⎜⎜

⎝

n i=0

1

n i=0

x_i

n i=0

x²_i

n i=0

xi

n i=0

x²_i

n i=0

x³_i

n i=0

x²_i

n i=0

x³_i

n i=0

x⁴_i

⎞

⎟⎟

⎠

⎛

⎜⎜

⎜⎝ c^∗₁ c^∗₂ c^∗₃

⎞

⎟⎟

⎟⎠=

⎛

⎜⎜

⎝

n i=0

f_i

n i=0

fixi

n i=0

fix²_i

⎞

⎟⎟

⎠ .

2. ϕ^∗(x) = −0.2835x²+ 1.2359x− 0.0130.

Shrnut´ı lekce

Uk´azali jsme z´akladn´ı postupy pro aproximaci funkc´ı (dat) pomoc´ı interpolace a metody nejmenˇs´ıch ˇctverc˚u.