5. Interpolace a aproximace funkcí

5. INTERPOLACE A APROXIMACE FUNKCÍ

5.

Numerické metody

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áme jej´ı funkˇcn´ı hodnoty fi = ´ lohy. f (xi ) v uzlech xi pro i = 0, . . . , n. Budeme rozliˇsovat dvˇe u Interpolaˇcn´ı u ´ loha: Hledáme funkci ϕ, pro niˇz je ϕ(xi ) = fi ,

i = 0, . . . , n.

(5.0.1)

Aproximace metodou nejmenˇs´ıch ˇctverc˚ u: Hledáme funkci ϕ, pro niˇz je ϕ(xi ) ≈ fi ,

i = 0, . . . , n,

(5.0.2)

kde pˇribliˇzná rovnost ≈” je urˇcena tak, aby souˇcet druchých mocnin odchylek ” ymi hodnotami ϕ(xi ) byl mimezi pˇredepsan´ ymi hodnotami fi a pˇredpokládan´ nimáln´ı.

Jestliˇze tyto u ´ lohy znázorn´ıme graficky, bude ˇreˇsen´ı interpolaˇcn´ı u ´ lohy procházet pˇres body (xi , fi ), i = 0, . . . , n, kdeˇzto ˇreˇsen´ı aproximaˇcn´ı u ´ lohy bude (obecnˇe) procházet jejich bl´ızk´ ym okol´ım. Formulace obou u ´ loh je zat´ım pˇr´ıliˇs obecná, protoˇze jsme neˇrekli jakého typu má b´ yt funkce ϕ. Ukáˇzeme tˇri volby: polynom, splajn (spline-funkce) a lineárn´ı kombinace obecn´ ych funkc´ı. Polynom je jednoduch´ y z hlediska matematick´ ych operac´ı (snadno se derivuje, integruje atp.), jeho graf vˇsak ˇcasto osciluje. Lepˇs´ı tvary grafu maj´ı splajny. Kombinace obecn´ ych 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.).

- 92 -


Numerické metody

5.1.

Interpolaˇ cn´ı polynom

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

Pˇredpokl´ adan´ e znalosti ˇ sen´ı soustav lineárn´ıch rovnic. Vˇeta o stˇredn´ı hodnotˇe difePolynomy. Reˇ renciáln´ıho poˇctu.

V´ yklad Funkci ϕ v u ´ loze (5.0.1) budeme hledat jako interpolaˇcn´ı polynom stupnˇe nejv´ yˇse n, tj. poloˇz´ıme ϕ = pn , kde pn (x) = a0 + a1 x + a2 x2 + . . . + an xn .

(5.1.1)

Zaˇcneme pˇr´ıkladem. Pˇ r´ıklad 5.1.1.

Jsou dány uzly x0 = −2, x1 = −1, x2 = 1, x3 = 2 a funkˇcn´ı

hodnoty f0 = 10, f1 = 4, f2 = 6, f3 = 3. Urˇcete interpolaˇcn´ı polynom p3 . ˇ sen´ı: Reˇ

Hledan´ y polynom má obecn´ y tvar p3 (x) = a0 + a1 x + a2 x2 + a3 x3 .

Koeficienty a0 , a1 , a2 , a3 urˇc´ıme tak, aby platilo (5.0.1). Kaˇzdá interpolaˇcn´ı rovnost urˇcuje jednu rovnici: p3 (−2) = 10 ⇒ a0 − 2a1 + 4a2 − 8a3 = 10, p3 (−1) = 4

⇒ a0 −

a1

+

a2

−

a3

=

4,

p3 (1)

= 6

⇒ a0 +

a1

+

a2

+

a3

=

6,

p3 (2)

= 3

⇒ a0 + 2a1 + 4a2 + 8a3 =

3.

- 93 -


Numerické metody

Dostali jsme soustavu lineárn´ıch rovnic ⎞ ⎞⎛ ⎞ ⎛ ⎛ a0 10 1 −2 4 −8 ⎟ ⎟⎜ ⎟ ⎜ ⎜ ⎟ ⎟⎜ ⎟ ⎜ ⎜ ⎜ 1 −1 1 −1 ⎟ ⎜ a1 ⎟ ⎜ 4 ⎟ ⎟ ⎟⎜ ⎟ ⎜ ⎜ ⎟, ⎟⎜ ⎟=⎜ ⎜ ⎜ ⎜ ⎟ ⎟ ⎜ 1 1 1 1 ⎟ ⎜ a2 ⎟ ⎜ 6 ⎟ ⎟ ⎜ ⎠ ⎠⎝ ⎠ ⎝ ⎝ 3 1 2 4 8 a3 jej´ımˇz ˇreˇsen´ım (na tˇri desetinná m´ısta) jsou koeficienty a0 = 4.500, a1 = 1.917, a2 = 0.500 a a3 = −0.917. Interpolaˇcn´ı polynom má tvar p3 (x) = 4.500 + 1.917x + 0.500x2 − 0.917x3 . Jeho graf je na obrázku 5.1.1. 12 10 8 6 4 2 −2

−1

0

1

2

Obrázek 5.1.1: Graf interpolaˇcn´ıho polynomu p3 . Rozborem postupu z pˇr´ıkaldu dokáˇzeme následuj´ıc´ı vˇetu. Vˇ eta 5.1.1. Necht’ jsou dány vzájemnˇe r˚ uzné uzly xi a funkˇcn´ı hodnoty fi , i = 0, . . . , n. Existuje právˇe jeden interpolaˇcn´ı polynom stupnˇe nejvýˇse n.

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

- 94 -


Numerické metody

(5.0.1) dostaneme soustavu lineárn´ıch rovnic pn (xi ) = a0 + a1 xi + a2 x2i + . . . + an xni = fi ,

i = 0, . . . , n,

kterou lze zapsat maticovˇe jako ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝

1 x0 x20 . . . xn0 1 x1

x21

...

xn1

1 x2 x22 . . . xn2 . .. .. .. . . . .. . . .

⎞⎛ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎠⎝

1 xn x2n . . . xnn

a0 a1 a2 .. . an

⎞

⎛

⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟=⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎠ ⎝

f0 f1 f2 .. .

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎟ ⎟ ⎠

fn

Matice této soustavy má nenulov´ y determinant (Vandermod˚ uv determinant). Odtud plyne existence jediného ˇreˇsen´ı soustavy lineárn´ıch rovnic a také inter2

polaˇcn´ıho polynomu.

5.1.1.

Lagrange˚ uv tvar interpolaˇ cn´ıho polynomu

Ukáˇzeme postup, pˇri nˇemˇz se obejdeme bez ˇreˇsen´ı soustavy lineárn´ıch rovnic. Interpolaˇcn´ı polynom budeme hledat ve tvaru pn (x) = f0 ϕ0 (x) + f1 ϕ1 (x) + . . . + fn ϕn (x).

(5.1.2)

Rovnosti pn (xi ) = fi , i = 0, 1, . . . , n budou splnˇeny, jestliˇze bude platit ⎧ ⎪ ⎨ 1 pro i = j, ϕi (xj ) = ⎪ ⎩ 0 pro i = j. Z vˇety 5.1.1. v´ıme, ˇze interpolaˇcn´ı polonom je stupnˇe nejvýˇse n, takˇze také vˇsechny funkce ϕi mus´ı b´ yt polynomy stupnˇe nejvýˇse n. Uvedeným poˇzadavk˚ um vyhovuje následuj´ıc´ı definice: ϕi (x) =

(x − x0 ) . . . (x − xi−1 )(x − xi+1 ) . . . (x − xn ) (xi − x0 ) . . . (xi − xi−1 )(xi − xi+1 ) . . . (xi − xn )

- 95 -

(5.1.3)


Numerické metody

ˇ pro i = 0, 1, . . . , n. Citatel je totiˇz polynom, kter´ y nab´ yvá nulov´ ych hodnot yvá nenulové hodnoty, která je ve vˇsech uzlech kromˇe xi . V uzlu xi pak nab´ obsaˇzena ve jmenovateli zlomku, takˇze plat´ı ϕi (xi ) = 1. aze interpolaˇcn´ı u ´ lohy Polynom˚ um ϕi , i = 0, 1, . . . , n se ˇr´ıká Lagrangeova b´ a vzorec (5.1.2) se naz´ yvá Lagrange˚ uv tvar inteprolaˇcn´ıho polynomu.

Pˇ r´ıklad 5.1.2.

Mˇejme dány uzly x0 = −2, x1 = −1, x2 = 1, x3 = 2 a funkˇcn´ı

uv tvar interpolaˇcn´ıho hodnoty f0 = 10, f1 = 4, f2 = 6, f3 = 3. Napiˇste Lagrange˚ polynomu. ˇ sen´ı: Reˇ

Nejdˇr´ıve sestav´ıme Lagrangeovu bázi. Podle (5.1.3) je

ϕ0 (x) =

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

ϕ1 (x) =

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

ϕ2 (x) =

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

ϕ3 (x) =

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

Dosazen´ım do (5.1.2) dostaneme výsledek 2 5 p3 (x) = − (x + 1)(x − 1)(x − 2) + (x + 2)(x − 1)(x − 2) − 6 3 1 −(x + 2)(x + 1)(x − 2) + (x + 2)(x + 1)(x − 1). 4 Pozn´ amka ´ Interpolaˇcn´ı polynom je podle vˇety 5.1.1. urˇcen jednoznaˇcnˇe. Upravou Lagrangeova tvaru proto mus´ıme nutnˇe doj´ıt k polynomu, který jsem vypoˇc´ıtali v pˇr´ıkladu 5.1.1. (ovˇeˇrte).

- 96 -


Numerické metody 5.1.2.

Newton˚ uv tvar interpolaˇ cn´ıho polynomu

Uvaˇzujme zápis polynomu ve tvaru: pn (x) = a0 +a1 (x−x0 )+a2 (x−x0 )(x−x1 )+. . .+an (x−x0 ) . . . (x−xn−1 ). (5.1.4) Jestliˇze dosad´ıme do interpolaˇcn´ıch rovnost´ı pn (xi ) = fi , i = 0, 1, . . . , n, dostaneme soustavu lineárn´ıch rovnic s doln´ı troj´ uheln´ıkovou matic´ı: ⎞⎛ ⎞ ⎛ ⎛ a 1 0 0 ... 0 ⎟⎜ 0 ⎟ ⎜ ⎜ ⎟⎜ ⎟ ⎜ ⎜ ⎜ ⎟ ⎜ ⎜ 1 (x1 − x0 ) 0 ... 0 ⎟ ⎟ ⎜ a1 ⎟ ⎜ ⎜ ⎟⎜ ⎟=⎜ ⎜ ⎜ 1 (x − x ) (x − x )(x − x ) . . . 0 ⎟ ⎜ a ⎟ ⎜ ⎟⎜ 2 ⎟ ⎜ ⎜ 2 0 2 0 2 1 ⎠⎝ ⎠ ⎝ ⎝ . .. .. .. .. .. . . . . . . .

⎞ f0

⎟ ⎟ f1 ⎟ ⎟ ⎟, f2 ⎟ ⎟ ⎠ .. .

(5.1.5)

Odud m˚ uˇzeme postupnˇe vyjádˇrit koeficienty ak : a0

a2

=

=

f0 ,

a1 =

f1 − a0 f1 − f0 = , x1 − x0 x1 − x0

f1 − f0 f2 − f1 − f2 − a1 (x2 − x0 ) − a0 x − x1 x1 − x0 , = 2 (x2 − x0 )(x2 − x1 ) x2 − x0

atd.. Výrazy na prav´ ych 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ány vzájemnˇe r˚ uzné uzly xi a funkˇcn´ı hodnoty fi , i = 0, . . . , n. Pomˇerné diference k-tého ˇra´du f [xi+k , . . . , xi ],

i = 0, 1, . . . , n − k definujeme

rekurentnˇe: • pro k = 0 : • pro k = 1 :

f [xi ] = fi ; f [xi+1 , xi ] =

• pro k ≤ n : f [xi+k , . . . , , xi ] =

fi+1 − fi ; xi+1 − xi f [xi+k , . . . , xi+1 ] − f [xi+k−1 , . . . , xi ] . xi+k − xi

- 97 -


Numerické metody

Porovnán´ım pomˇern´ ych diferenc´ı s koeficienty ak vid´ıme, ˇze ak = f [xk , . . . , x0 ],

k = 0, 1, . . . , n.

Dosazen´ım do (5.1.4) dostaneme Newton˚ uv tvar interpolaˇcn´ıho polynomu: pn (x) = f0 + f [x1 , x0 ](x − x0 ) + . . . + f [xn , . . . , x0 ](x − x0 ) . . . (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.


uv tvar interpolaˇcn´ıho hodnoty f0 = 10, f1 = 4, f2 = 6, f3 = 3. Napiˇste Newton˚ polynomu. ˇ sen´ı: Reˇ

Potˇrebujeme vypoˇc´ıtat pomˇerné diference: f [x1 , x0 ], f [x2 , x1 , x0 ], f [x3 , x2 , x1 , x0 ].

Podle definice je f [x1 , x0 ] =

4 − 10 = −6, −1 + 2

f [x2 , x1 ] − f [x1 , x0 ] f [x2 , x1 , x0 ] = = x2 − x0

6−4 1+1

+6 7 = , 1+2 3

f [x3 , x2 , x1 ] − f [x2 , x1 , x0 ] f [x3 , x2 , x1 , x0 ] = = x3 − x0

3−6 − 6−4 2−1 1+1

− 2+2

2+1

7 3

= −

11 . 12

Dosazen´ım do (5.1.6) dostaneme výsledek 11 7 p3 (x) = 10 − 6(x + 2) + (x + 2)(x + 1) − (x + 2)(x + 1)(x − 1). 3 12 Pˇrehlednˇe m˚ uˇzeme výpoˇcet pomˇern´ ych 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 ˇra´d˚ u. Pro napsán´ı interpolaˇcn´ıho polynomu potˇrebujeme z této tabulky hodnoty diferenc´ı z prvn´ıho ˇra´dku.

- 98 -


Numerické metody

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

5.1.3.

xi

fi

f [xi+1 , xi ] f [xi+2 , xi+1 , xi ] f [x3 , x2 , x1 , x0 ] −6

7 3

4

1

− 43

1

6

−3

2

3

0

−2 10

1

−1

2 3

− 11 12

Interpolaˇ cn´ı chyba

Pˇredpokládejme, ˇze hodnoty fi jsou funkˇcn´ımi hodnotami funkce f v uzlech xi , tj. fi = f (xi ). Bude nás zaj´ımat interpolaˇcn´ı chyba f (x) − pn (x). uˇze být velká. V uzlech xi je interpolaˇcn´ı chyba nulová, ale mimo uzly m˚ Vˇ eta 5.1.2. uzné a leˇz´ı na intervalu a, b. Necht’ Necht’ uzly xi , i = 0, 1, . . . , n, jsou vzájemnˇe r˚ funkce f má na tomto intervalu n + 1 spojit´ ych derivac´ı. Pak pro kaˇzdé x ∈ a, b existuje ξ = ξ(x) v (a, b) tak, ˇze plat´ı f (x) − pn (x) =

f (n+1) (ξ) πn+1 (x), (n + 1)!

(5.1.7)

kde πn+1 (x) = (x − x0 ) . . . (x − xn ). D˚ ukaz: Pro x = xi je rovnost (5.1.7) splnˇena, protoˇze obˇe jej´ı strany jsou nulové. Pro pevnˇe zvolené x = xi definujme funkci g(t) = f (t) − pn (t) −

πn+1 (t) (f (x) − pn (x)) , πn+1 (x)

(5.1.8)

kde t je promˇenná a x je parametr. Funkce g má zˇrejmˇe n+2 koˇren˚ u, kter´ ymi jsou body x0 , . . ., xn a x. Kaˇzdá derivace funkce g má 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

- 99 -

5. INTERPOLACE A APROXIMACE FUNKCÍ (n+1)

v´ yraz (5.1.8) (podle t) a pouˇzijeme pˇritom pn

Numerické metody (n+1)

(t) = 0 a πn+1 (t) = (n + 1)!,

dostaneme 0 = g (n+1) (ξ) = f (n+1) (ξ) −

(n + 1)! (f (x) − pn (x)) . πn+1 (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´ y vliv tvar polynomu πn+1 , 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 + x2

a graf interpolaˇcn´ıho polynomu odpov´ıdaj´ıc´ıho uzl˚ um xi = −5+i, i = 0, 1, . . . , 10. Výsledek porovnáme s grafem polynomu π11 (x) = (x + 5)(x + 4) . . . (x − 5). ˇ sen´ı: Obrázek 5.1.2.a ukazuje graf polynomu π11 . Z jeho pr˚ ubˇehu lze usouReˇ dit, ˇze nejvˇetˇs´ı interpolaˇcn´ı chyby budou pobl´ıˇz krajn´ıch uzl˚ u x0 = −5 a x10 = 5. Na obrázku 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ávˇe na kraj´ıch intervalu −5, 5. Poznamenejme 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ázka 2. Jakého stupnˇe je interpolaˇcn´ı polynom? Otázka 3. Jak se chová interpolaˇcn´ı chyba?

´ Ulohy k samostatn´ emu ˇreˇsen´ı 1. Pro uzly x0 = −1, x1 = 0, x2 = 2, x3 = 3, x4 = 5 a funkˇcn´ı hodnoty f0 = −2,

- 100 -


Numerické metody 5

5

x 10

2 1.5 1

0 0.5 0 −5 −5

0

−0.5 −5

5

a

0

5

b

Obrázek 5.1.2: a) Graf π11 ; b) Grafy f (neosciluj´ıc´ı) a p10 (osciluj´ıc´ı). f1 = 1, f2 = 0, f3 = 2, f4 = −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´ ysledky u ´loh k samostatn´ emu ˇreˇsen´ı 3 4 x + 1. p4 (x) = − 20

11 3 x 10

−

109 2 x 60

−

1 x 15

+ 1.

1 1 2. Lagrange˚ uv tvar: p4 (x) = − 36 x(x − 2)(x − 3)(x − 5) − 30 (x + 1)(x − 2)(x − 3)

(x − 5) −

1 (x 12

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

1 (x 180

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

3 Newton˚ uv tvar: p4 (x) = −2 + 3(x + 1) − 35 (x + 1)x + 12 (x + 1)x(x − 2) − 20 (x + 1) 30

x(x − 2)(x − 3).

- 101 -


5.2.

Numerické metody

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´ı u ´ lohy pouˇz´ıt funkci, která bude poˇcástech polynomem n´ızkého stupnˇe, jej´ıˇz jednotlivé ˇcásti budou na sebe navazovat dostateˇcnˇe hladce. Takov´ ym funkc´ım se ˇr´ıká splajn (z angl. spline”). Ukáˇzeme dva nejˇcastˇeji pouˇz´ıvané splajny: lineárn´ı a kubický. ” Pˇredpokl´ adan´ e znalosti ˇ sen´ı soustav lineárn´ıch rovnic. Interpolaˇcn´ı polynom. Spojitost derivace. Reˇ

V´ yklad Abychom se vyhnuli komplikac´ım pˇri popisu, budeme pˇredpokládat, ˇze uzly interpolace tvoˇr´ı rostouc´ı posloupnost, tzn. x0 < x1 < . . . < xn . Vzdálenost dvou sousedn´ıch uzl˚ u oznaˇc´ıme hi , tj. hi = xi − xi−1 , 5.2.1.

i = 1, . . . , n.

Line´ arn´ı splajn

Definice 5.2.1. Lineárn´ım splajnem nazýváme funkci s1 , která je spojitá na intervalu x0 , xn a na kaˇzdém podintervalu xi−1 , xi , i = 1, . . . , n, je polynomem prvn´ıho stupnˇe. Lineárn´ı interpolaˇcn´ı splajn je ˇreˇsen´ım u ´ lohy (5.0.1), tzn. ˇze pro nˇej plat´ı uˇzeme ho zapsat po ˇcástech pro i = 1, . . . , n: s1 (xi ) = fi , i = 0, . . . , n. M˚ s1 (x) = fi−1 (1 − t) + fi t,

t = (x − xi−1 )/hi ,

Grafem lineárn´ıho splajnu je lomen´ a ˇcára.

- 102 -

x ∈ xi−1 , xi .

(5.2.1)


Numerické metody

Pˇ r´ıklad 5.2.1.


hodnoty f0 = 10, f1 = 4, f2 = 6, f3 = 3. Napiˇste lineárn´ı interpolaˇcn´ı splajn. ˇ sen´ı: Zap´ıˇseme jej pomoc´ı pˇredpisu (5.2.1): Reˇ ⎧ ⎪ 10ϕ1 (t) + 4ϕ2 (t), t = x + 2 ⎪ ⎪ ⎪ ⎨ s1 (x) = 4ϕ1 (t) + 6ϕ2 (t), t = (x + 1)/2 ⎪ ⎪ ⎪ ⎪ ⎩ 6ϕ1 (t) + 3ϕ2 (t), t = x − 1

pro x ∈ −2, −1, pro x ∈ −1, 1, pro x ∈ 1, 2,

kde ϕ1 (t) = 1 − t, ϕ2 (t) = t. Graf je znázornˇen na obrázku 5.2.1. 5.2.2.

Kubick´ y splajn

Definice 5.2.2. Kubickým splajnem naz´ yváme funkci s3 , která má na intervalu x0 , xn dvˇe spojité derivace a na kaˇzdém podintervalu xi−1 , xi , i = 1, . . . , n, je polynomem tˇret´ıho stupnˇe. Kubický interpolaˇcn´ı splajn, je ˇreˇsen´ı interpolaˇcn´ı u ´ lohy (5.0.1). Jeho konstrukce 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: s3 (x) = fi−1 (1 − 3t2 + 2t3 ) + fi (3t2 − 2t3 ) +mi−1 hi (t − 2t2 + t3 ) + mi hi (−t2 + t3 ),

(5.2.2)

kde t = (x − xi−1 )/hi a x ∈ xi−1 , xi . Tento pˇredpis je navrˇzen tak, aby parameyznam funkˇcn´ıch hodnot a hodnot prvn´ı derivace try fi−1 , fi a mi−1 , mi mˇely v´ v krajn´ıch bodech intervalu xi−1 , xi , tj. plat´ı s3 (xi−1 ) = fi−1 , s3 (xi−1 ) = mi−1 ,

- 103 -

s3 (xi ) = fi , s3 (xi ) = mi .

(5.2.3) (5.2.4)


Numerické metody

O splnˇen´ı rovnost´ı (5.2.3) a (5.2.4) se m˚ uˇzeme pˇresvˇedˇcit dosazen´ım xi−1 a xi do (5.2.2) a do prvn´ı derivace s3 , kterou vyjádˇr´ıme z (5.2.2) podle pravidla o derivován´ı sloˇzené funkce: s3 (x) = fi−1 (−6t + 6t2 )/hi + fi (6t − 6t2 )/hi +mi−1 (1 − 4t + 3t2 ) + mi (−2t + 3t2 ).

(5.2.5)

Pˇredpis (5.2.2) zaruˇcuje spojitost prvn´ı derivace s3 na celém intervalu x0 , xn pro libovolné hodnoty mi . Spojitost druhé derivace vynut´ıme speciáln´ı volbou mi . Budeme poˇzadovat lim s3 (x) = lim s3 (x)

x→xi −

(5.2.6)

x→xi +

yraz pro druhou derivaci ve vnitˇrn´ıch uzlech xi , i = 1, . . . , n − 1. Potˇrebný v´ vypoˇcteme z (5.2.5) opˇet podle pravidla o derivován´ı sloˇzené funkce: s3 (x) = fi−1 (−6 + 12t)/h2i + fi (6 − 12t)/h2i +mi−1 (−4 + 6t)/hi + mi (−2 + 6t)/hi .

(5.2.7)

Levou stranu v (5.2.6) vyjádˇr´ıme z (5.2.7) pro t = 1: lim s3 (x) = 6fi−1 /h2i − 6fi /h2i + 2mi−1 /hi + 4mi /hi .

x→xi −

(5.2.8)

Pravou stranu v (5.2.6) vyjádˇr´ıme z (5.2.7) pro t = 0, kdyˇz souˇcasnˇe posuneme indexován´ı: lim s (x) = −6fi /h2i+1 + 6fi+1 /h2i+1 − 4mi /hi+1 − 2mi+1 /hi+1 .

x→xi +

(5.2.9)

Dosad´ıme-li (5.2.8) a (5.2.9) do (5.2.6), dostaneme po jednoduché u ´ pravˇe hi+1 mi−1 + 2(hi+1 + hi )mi + hi mi+1 =

! hi+1 hi hi hi+1 fi−1 + − fi+1 , fi + 3 − hi hi hi+1 hi+1

- 104 -

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


Numerické metody

Tyto rovnosti tvoˇr´ı soustavu n − 1 rovnic pro n + 1 neznám´ ych mi , i = 0, 1, . . . .n. Abychom dostali jediné ˇreˇsen´ı, urˇc´ıme m0 a mn napˇr´ıklad jako pˇribliˇzné derivace: m0 = Pˇ r´ıklad 5.2.2.

f1 − f0 , h1

mn =

fn − fn−1 . hn

(5.2.11)


hodnoty f0 = 10, f1 = 4, f2 = 6, f3 = 3. Napiˇste kubick´ y interpolaˇcn´ı splajn. ˇ sen´ı: Reˇ

Nejdˇr´ıve vypoˇc´ıtáme parametry mi , i = 0, 1, 2, 3. Podle (5.2.11) je m0 =

4 − 10 = −6, −1 + 2

m3 =

3−6 = −3. 2−1

Soustava (5.2.10) má dvˇe rovnice: 2(h2 + h1 )m1 + h1 m2 h3 m1 + 2(h3 + h2 )m2 které m˚ uˇzeme psát jako ⎛

! h1 f1 + f2 − h2 m0 , h2

! h3 h3 h2 h2 = 3 − f1 + − f2 + f3 − h2 m3 , h2 h2 h3 h3 h2 = 3 − f0 + h1

⎞⎛

⎞

h2 h1 − h1 h2

⎛

⎞

⎜ 6 1 ⎟ ⎜ m1 ⎟ ⎜ −21 ⎟ ⎝ ⎠⎝ ⎠=⎝ ⎠. m2 1 6 −9

, m2 = − 66 . Výsledn´ y splajn zap´ıˇseme podle Vyˇreˇsen´ım dostaneme m1 = − 234 70 70 (5.2.2) po ˇcástech: ⎧ ⎪ 10ϕ1 (t) + 4ϕ2 (t) − 6ϕ3 (t) − 117 ϕ (t), ⎪ 35 4 ⎪ ⎪ ⎨ s3 (x) = 4ϕ1 (t) + 6ϕ2 (t) − 234 ϕ (t) − 66 ϕ (t), 35 3 35 4 ⎪ ⎪ ⎪ ⎪ ⎩ ϕ (t) − 3ϕ4 (t), 6ϕ1 (t) + 3ϕ2 (t) − 33 35 3

t = x + 2 pro x ∈ −2, −1, t = (x + 1)/2 pro x ∈ −1, 1, t = x − 1 pro x ∈ 1, 2,

kde ϕ1 (t) = 1 − 3t2 + 2t3 , ϕ2 (t) = 3t2 − 2t3 , ϕ3 (t) = t − 2t2 + t3 , ϕ4 (t) = −t2 + t3 . Graf je znázornˇen na obrázku 5.2.1. Pˇ r´ıklad 5.2.3.

(Rungeho pˇ r´ıklad, pokraˇ cov´ an´ı) Nakresl´ıme graf inter-

polaˇcn´ıho kubického splajnu pro funkci f a uzly xi z pˇr´ıkladu 5.1.4. a porovnáme ho s grafem interpolaˇcn´ıho polynomu.

- 105 -


Numerické metody

10 8 6

s1 4 2

s3 −2

−1

0

1

2

Obrázek 5.2.1: Graf lineárn´ıho (s1 ) a kubického interpolaˇcn´ıho (s3 ) splajnu. ˇ sen´ı: Reˇ

Na obrázku 5.2.2 vid´ıme, ˇze splajn s3 neosciluje a je proto mno-

hem lepˇs´ı aproximac´ı interpolované funkce f neˇz interpolaˇcn´ı polynom p10 , viz obrázek 5.1.2.b. 1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

−5

0

0

5

a

−5

0

5

b

Obrázek 5.2.2: a) Funkce f ; b) Kubick´ y interpolaˇcn´ı splajn s3 .

Kontroln´ı ot´ azky Otázka 1. Co je to splajn? Jak se definuje a poˇc´ıtá splajn lineárn´ı a kubick´ y? Otázka 2. Jak se chovaj´ı pˇri interpolaci splajny v porovnán´ı s polynomy?

- 106 -


Numerické metody

´ Ulohy k samostatn´ emu ˇreˇsen´ı 1. Pro uzly x0 = −1, x1 = 0, x2 = 2, x3 = 3, x4 = 5 a funkˇcn´ı hodnoty f0 = −2, f1 = 1, f2 = 0, f3 = 2, f4 = −1 sestavte lineárn´ı interpolaˇcn´ı splajn. 2. Pro stejná data sestavte kubick´ y interpolaˇcn´ı splajn. V´ ysledky u ´loh k samostatn´ emu ˇreˇsen´ı 1.

s1 (x) =

⎧ ⎪ −2ϕ1 (t) + ϕ2 (t), t = x + 1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ϕ1 (t), t = x/2 ⎪ ⎪ ⎪ 2ϕ2 (t), ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 2ϕ (t) − 1ϕ (t), 1 2

pro x ∈ −1, 0, pro x ∈ 0, 2,

t=x−2

pro x ∈ 2, 3,

t = (x − 3)/2

pro x ∈ 3, 5,

kde ϕ1 (t) = 1 − t, ϕ2 (t) = t. 2. Krajn´ı parametry jsou m0 = 3, m4 = − 32 , ostatn´ı dostaneme ze soustavy ⎛

6 1 0

⎜ ⎜ ⎜ 1 6 2 ⎜ ⎝ 0 2 6 takˇze m1 =

s3 (x) =

97 , 62

m2 =

69 , 62

m3 =

35 31

⎞⎛

m1

⎞

⎛

21 2

⎟⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ m ⎟ = ⎜ 21 ⎟⎜ 2 ⎟ ⎜ 2 ⎠⎝ ⎠ ⎝ m3 9

⎞ ⎟ ⎟ ⎟, ⎟ ⎠

a koneˇcnˇe

⎧ ⎪ −2ϕ1 (t) + ϕ2 (t) + 3ϕ3 (t) + ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ϕ1 (t) + 97 ϕ (t) + 69 ϕ (t), 31 3 31 4

97 ϕ (t), 62 4

⎪ ⎪ ⎪ ϕ (t) + 35 3ϕ4 (t), 2ϕ2 (t) + 69 ⎪ 62 3 31 ⎪ ⎪ ⎪ ⎪ ⎩ 2ϕ (t) − ϕ (t) + 70 ϕ (t) − 3ϕ (t), 1 2 4 31 3

t = x + 1 pro x ∈ −1, 0, t = x/2 pro x ∈ 0, 2, t = x − 2 pro x ∈ 2, 3, t = (x − 3)/2 pro x ∈ 3, 5,

kde ϕ1 (t) = 1 − 3t2 + 2t3 , ϕ2 (t) = 3t2 − 2t3 , ϕ3 (t) = t − 2t2 + t3 , ϕ4 (t) = −t2 + t3 .

- 107 -


5.3.

Numerické metody

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. Interpolace je nepouˇzitelná, jestliˇze je poˇzadován jist´ y 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.

Pˇredpokl´ adan´ e znalosti ˇ sen´ı Lineárn´ı závislost a nezávislost. Urˇcen´ı minima funkce pomoc´ı derivace. Reˇ soustav lineárn´ıch rovnic. V´ yklad Zaˇcneme pˇr´ıkladem. Pˇ r´ıklad 5.3.1. Mˇejme dány uzly x0 = −2, x1 = −1, x2 = 1, x3 = 2 a funkˇcn´ı hodnoty f0 = 10, f1 = 4, f2 = 6, f3 = 3. Najdˇete pˇr´ımku ϕ(x) = c1 + c2 x,

(5.3.1)

která je bl´ızko” pˇredepsan´ ym hodnotám. ” ˇ sen´ı: Reˇ

Nejdˇr´ıve se mus´ıme rozhodnout jak chápat slovo bl´ızko”. Uˇz jsme ” to vlastnˇe ˇrekli, kdyˇz jsme popisovali smysl pˇribliˇzných rovnost´ı v aproximaˇcn´ı u ´ loze (5.0.2). Pˇr´ımku ϕ urˇc´ıme tak, aby minimalizovala souˇcet druhých mocnin ´ lohu na miodchylek 3i=0 (ϕ(xi )−fi )2 . Jestliˇze sem dosad´ıme (5.3.1), dostaneme u 3 nimalizaci funkce dvou promˇenných Ψ(c1 , c2 ) = i=0 (c1 + c2 xi − fi )2 . Minimum

- 108 -


Numerické metody c∗1 , c∗2 vyhovuje rovnic´ım

∂Ψ ∗ ∗ (c , c ) = 0, ∂c1 1 2

∂Ψ ∗ ∗ (c , c ) = 0. ∂c2 1 2

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

3

(c∗1 + c∗2 xi − fi ) = 0,

2

i=0

3

(c∗1 + c∗2 xi − fi )xi = 0,

i=0

coˇz je soustava lineárn´ıch rovnic ⎛

3

3

i=0

i=0

1 xi ⎜ ⎜ i=0 i=0 ⎜ 3 3 ⎜ ⎝ xi x2i

⎞

⎛

⎞

⎛

⎞

3

fi ⎜ ⎟ ⎜ ⎟ c∗1 i=0 ⎟⎝ ⎠=⎜ 3 ⎜ ⎟ ⎝ ⎠ c∗2 fi xi

⎞ ⎞ ⎛ ⎞⎛ ⎛ ⎟ ∗ ⎟ 23 c 4 0 ⎟ , tj. ⎝ ⎠. ⎠⎝ 1 ⎠ = ⎝ ⎟ ⎠ −12 c∗2 0 10

i=0

Tato soustava má jediné ˇreˇsen´ı c∗1 =

23 , 4

c∗2 = − 65 , takˇze hledaná pˇr´ımka ϕ∗ = ϕ

je urˇcena pˇredpisem ϕ∗ (x) =

23 6 − x. 4 5

(5.3.2) 2

Jej´ı graf je znázornˇen na obrázku 5.3.1.

10 8 6 4 2 −2

−1

0

1

2

Obrázek 5.3.1: Aproximace metodou nejmenˇs´ıch ˇctverc˚ u; pˇr´ımka (5.3.2) plnˇe; funkce (5.3.8) ˇcárkovanˇe.

- 109 -


Numerické metody

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) = c1 ϕ1 (x) + . . . + cm ϕm (x) =

m

cj ϕj (x),

(5.3.3)

j=1

kde koeficienty c1 , . . . , cm jsou libovolná ˇc´ısla. Funkci ϕ∗ , pro niˇz plat´ı n

(ϕ∗ (xi ) − fi )2 ≤

i=0

n

(ϕ(xi ) − fi )2

∀ϕ

(5.3.4)

i=0

naz´ yváme aproximac´ı podle metody nejmenˇs´ıch ˇctverc˚ u. Jej´ı koeficienty c∗1 , . . . , c∗m urˇc´ıme jako minimum funkce n m Ψ(c1 , . . . , cm ) = ( cj ϕj (xi ) − fi )2 ,

(5.3.5)

i=0 j=1

které vyhovuje rovnic´ım ∂Ψ ∗ (c , . . . , c∗m ) = 0, ∂ck 1

k = 1, . . . , m.

(5.3.6)

Vyjádˇr´ıme-li parciáln´ı derivace n m ∂Ψ =2 ( cj ϕj (xi ) − fi )ϕk (xi ) ∂ck i=0 j=1

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

n m j=1

i=0

ϕj (xi )ϕk (xi )

c∗j

=

n

fi ϕk (xi ),

k = 1, . . . , m.

(5.3.7)

i=0

Soustava (5.3.7) se naz´ yvá soustava norm´ aln´ıch rovnic. Vˇ eta 5.3.1. Necht’ jsou dány vzájemnˇe r˚ uzné uzly xi a funkˇcn´ı hodnoty fi , i = 0, . . . , n. Necht’ je dán systém funkc´ı ϕj , j = 1, . . . , m, které jsou lineárnˇe nezávislé. Potom nuje (5.3.4) a jej´ı koeficienty c∗1 , . . . , c∗m jsou existuje jediná funkce ϕ∗ , která splˇ ˇreˇsen´ım soustavy normáln´ıch rovnic (5.3.7).

- 110 -


Numerické metody

D˚ ukaz: V bodˇe c∗1 , . . . , c∗m , kter´ y vyhovuje rovnic´ım (5.3.6), nab´ yvá funkce Ψ minima, jestliˇze matice druhých derivac´ı je symmetrická a pozitivnˇe definitn´ı (kladná). Druhé derivace jsou urˇceny vzorcem ∂2Ψ =2 ϕl (xi )ϕk (xi ), ∂ck ∂cl i=0 n

odkud je symetrie zˇrejmá na prvn´ı pohled (prohozen´ım index˚ u k a l se nic nezmˇen´ı). Necht’ d1 , . . . , dm jsou ˇc´ısla ne vˇsechny souˇcasnˇe nulová. Potom n " m m n #" # ∂2Ψ dk dl =2 dl ϕl (xi ) dk ϕk (xi ) = 2 ϕ(x ˜ i )2 > 0, ∂c k ∂cl i=0 l=1 i=0 k=1 l=1 k=1 m kde ϕ(x) ˜ = k=1 dk ϕk (x), takˇze matice druhých derivac´ı je pozitivnˇe definitn´ı. m m

Odtud také plyne, ˇze matice soustavy normáln´ıch rovnic je regulárn´ı, coˇz znamená, ˇze existuje jej´ı jediné ˇreˇsen´ı c∗1 , . . . , 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´ y lineárnˇe nezávisl´ y systém 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´ı ϕ1 (x) = e−x ,

ϕ2 (x) = sin x.

Aproximujte data z pˇr´ıkladu 5.3.1. ˇ sen´ı: Obecnˇe má soustava normáln´ıch rovnic tvar Reˇ ⎛ ⎞ ⎛ n n n ⎛ ⎞ −2xi −xi e e sin x fi e−xi i ⎟ ∗ ⎜ ⎜ c ⎟ ⎝ 1 ⎠ ⎜ i=0 ⎜ i=0 i=0 =⎜ ⎟ ⎜ n n n ⎝ ⎠ c∗ ⎝ −xi 2 2 e sin xi sin xi fi sin xi i=0

i=0

i=0

Po dosazen´ı dostaneme ⎛ ⎞⎛ ⎞ ⎞ ⎛ ∗ ⎜ 62.1409 −8.5736 ⎟ ⎝ c1 ⎠ ⎜ 87.3770 ⎟ =⎝ ⎝ ⎠ ⎠ ∗ c −8.5736 3.0698 −4.6821 2

- 111 -

⎞ ⎟ ⎟ ⎟. ⎠


Numerické metody

a odtud vypoˇc´ıtáme c∗1 = 1.9452, c∗2 = 3.9076, tj. ϕ∗ (x) = 1.9452e−x + 3.9076 sin x.

(5.3.8)

Graf je znázornˇen na obrázku 5.3.1. Kontroln´ı ot´ azky Otázka 1. Kdy je vhodné 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ázka 4. Co se stane, kdyˇz v (5.3.3) a (5.3.4) bude m = n + 1? ´ Ulohy k samostatn´ emu ˇreˇsen´ı 1. Napiˇste soustavu normáln´ıch lineárn´ıch rovnic pro systém funkc´ı ϕ1 (x) = 1, ϕ2 (x) = x, ϕ3 (x) = x2 . 2. Data x0 = −1, x1 = 0, x2 = 2, x3 = 3, x4 = 5 a f0 = −2, f1 = 1, f2 = 0, f3 = 2, f4 = −1 aproximujte metodou nejmenˇc´ıch ˇctverc˚ u pomoc´ı systému funkc´ı z pˇredchoz´ı u ´ lohy. V´ ysledky u ´loh k samostatn´ emu ˇreˇsen´ı 1. Normáln´ı soustava má tvar: ⎛ n ⎛ n ⎞ ⎞ n n 2 1 xi xi ⎟ ⎛ fi ⎟ ⎞ ⎜ ⎜ ⎜ i=0 ⎜ i=0 ⎟ c∗ ⎟ i=0 i=0 ⎜ ⎟⎜ 1 ⎟ ⎜ ⎟ n n n ⎜ ⎜ n ⎟ ⎟ 2 3 ⎟⎜ ∗ ⎟ ⎜ ⎜ ⎟. x x x f x = ⎟ ⎜ i i i c i i ⎜ ⎟⎝ 2 ⎠ ⎜ ⎟ ⎜ i=0 ⎜ i=0 ⎟ ⎟ i=0 i=0 ⎜ ⎜ ⎟ c∗ ⎟ n n n n 3 ⎝ ⎝ ⎠ ⎠ x2i x3i x4i fi x2i i=0

i=0

i=0

i=0

2. ϕ∗ (x) = −0.2835x2 + 1.2359x − 0.0130. Shrnut´ı lekce Ukázali jsme základn´ı postupy pro aproximaci funkc´ı (dat) pomoc´ı interpolace a metody nejmenˇs´ıch ˇctverc˚ u.

- 112 -

5. Interpolace a aproximace funkcí

Recommend Documents