Aproximace funkcí. Numerické metody. Ústav matematiky. 1. listopadu 2006

Aproximace funkc´ı

Numerické metody ˇ Doc. RNDr. Libor Cerm´ ak, CSc.

RNDr. Rudolf Hlaviˇcka, CSc.

´ Ustav matematiky Fakulta strojn´ıho inˇzen´ yrstv´ı Vysok´ e uˇ cen´ı technick´ e v Brnˇ e

1. listopadu 2006


Obsah

3

Aproximace funkc´ı Aproximace a interpolace Interpolace Interpolace polynomem Interpolaˇ cn´ı splajny Interpolace funkc´ı v´ıce promˇ enn´ ych

Metoda nejmenˇs´ıch ˇctverc˚ u Literatura

Aproximace a interpolace



Aproximovat funkci f (x) znamená nahradit ji funkc´ı ϕ(x), která je k f (x) v jistém smyslu bl´ızká. P´ıˇseme ϕ(x) ≈ f (x). Budeme se zabývat dvˇema základn´ımi typy aproximace, a to interpolac´ı a metodou nejmenˇs´ıch ˇctverc˚ u. Interpolace je taková aproximace, pˇri n´ıˇz ϕ(x) nabývá v zadaných bodech xi pˇredepsaných hodnot yi = f (xi ). Nˇekdy nav´ıc ˇzádáme, aby funkce ϕ a f mˇely v bodech xi také stejné derivace. Interpolaci je vˇenován odstavec 5. Metoda nejmenˇs´ıch ˇctverc˚ u je taková aproximace, pˇri n´ıˇz ϕ(x) prokládáme“ ” mezi zadanými body [xi , yi ] tak, aby vzdálenost“ funkc´ı f a ϕ byla ” v jistém smyslu minimáln´ı. Je pˇritom charakteristické, ˇze funkce ϕ body [xi , yi ] neprocház´ı. Metoda nejmenˇs´ıch ˇctverc˚ u je vyloˇzena v odstavci 50.



Aproximaci ϕ(x) pouˇzijeme k pˇribliˇznému výpoˇctu hodnot funkce f (x), tˇreba pˇri vykreslován´ı ϕ ≈ f . Je ˇzádouc´ı, aby výpoˇcet ϕ(x) byl jednoduchý“. Proto se ϕ ” ˇcasto hledá ve tvaru polynomu. Obecnˇe, ϕ(x) se pouˇz´ıvá k ˇreˇsen´ı u ´loh, v nichˇz vystupuje funkce f , kterou je u ´ˇcelné nebo dokonce nezbytné nahradit jej´ı vhodnou aproximac´ı ϕ. Jako pˇr´ıklad uved’me výpoˇcet derivace nebo urˇcitého integrálu: f 0 (x) nahrad´ıme pomoc´ı ϕ0 (x) Rb Rb a a f (x) dx nahrad´ıme pomoc´ı a ϕ(x) dx.


Obsah

3



Interpolace


Interpolace

Interpolaˇcn´ı funkci ϕ(x) vyb´ıráme z vhodné tˇr´ıdy funkc´ı. Omez´ıme se na dva nejbˇeˇznˇejˇs´ı pˇr´ıpady: a) ϕ(x) je polynom; b) ϕ(x) je po ˇcástech polynom, na kaˇzdém subintervalu obecnˇe jiný.


Interpolace

Interpolace polynomem

Pˇredpokládejme, ˇze jsou dány navzájem r˚ uzné body x0 , x1 , . . . , xn ,

xi 6= xj

pro i 6= j ,

ˇr´ıkáme jim také uzly interpolace, a v kaˇzdém z nich je pˇredepsána hodnota yi . Hledáme interpolaˇcn´ı polynom Pn (x) stupnˇe nejvýˇse n, který splˇ nuje interpolaˇcn´ı podm´ınky Pn (xi ) = yi , i = 0, 1, . . . , n . (3.1)


Interpolace

Existenci interpolaˇcn´ıho polynomu dokáˇzeme tak, ˇze ho zkonstruujeme. Lagrange˚ uv tvar interpolaˇ cn´ıho polynomu má vyjádˇren´ı Pn (x) = y0 `0 (x) + y1 `1 (x) + · · · + yn `n (x) =

n X

yi ì (x) ,

(3.2)

i=0

kde ì (x) jsou tzv. fundamentáln´ı polynomy definované pˇredpisem ì (x) =

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

(3.3)

Snadno nahlédneme, ˇze ì (xk ) =

1 0

pro k = i, pro k = 6 i,

takˇze interpolaˇcn´ı podm´ınky Pn (xk ) = splnˇeny.

Pn

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

i=0 yi ì (xk )

(3.4)

= yk , k = 0, 1, . . . , n, jsou


Interpolace

Interpolaˇcn´ı polynom je daty [xi , yi ], i = 0, 1, . . . , n, urˇcen jednoznaˇcnˇe. Skuteˇcnˇe, jsou-li P a Q interpolaˇcn´ı polynomy splˇ nuj´ıc´ı interpolaˇcn´ı podm´ınky P(xi ) = Q(xi ) = yi , pak polynom P − Q je roven nule v uzlech x0 , x1 , . . . , xn . Avˇsak polynom stupnˇe nejvýˇse n nem˚ uˇze m´ıt v´ıce neˇz n koˇren˚ u, pokud se nerovná identicky nule. V naˇsem pˇr´ıpadˇe proto nutnˇe P − Q = 0 a tedy P = Q. T´ım je jednoznaˇcnost interpolaˇcn´ıho polynomu dokázána.


Interpolace

Pˇr´ıklad 3.1. Urˇc´ıme interpolaˇcn´ı polynom pro data pˇredepsaná tabulkou xi

−1

1

2

3

yi

−6

−2

−3

2

Nejdˇr´ıve z´ıskáme fundamentáln´ı polynomy `0 (x) =

(x − 1)(x − 2)(x − 3) 1 = − (x 3 − 6x 2 + 11x − 6) , (−1 − 1)(−1 − 2)(−1 − 3) 24

`1 (x) =

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

`2 (x) =

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

`3 (x) =

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

a pak sestav´ıme interpolaˇcn´ı polynom P3 (x) = −6 · `0 (x) − 2 · `1 (x) − 3 · `2 (x) + 2 · `3 (x) = x 3 − 3x 2 + x − 1 .

2


Interpolace

Hlavn´ı pˇrednost´ı Lagrangeova tvaru interpolaˇcn´ıho polynomu je jeho elegantn´ı forma. Pouˇz´ıvá se proto zejména v teoretických u ´vahách. Pro praktické pouˇzit´ı vˇsak ideáln´ı nen´ı. Upozornˇeme na dva jeho hlavn´ı nedostatky. (a) Pˇridáme-li dalˇs´ı uzel xn+1 , mus´ıme pˇrepoˇc´ıtat vˇsechny fundamentáln´ı polynomy. (b) Poˇcet operac´ı potˇrebných k výpoˇctu hodnoty Pn (¯ x ) je pomˇernˇe znaˇcný, vyˇzaduje 2n2 + 2n násobic´ıch operac´ı a 2n2 + 3n operac´ı sˇc´ıtac´ıch. Oba tyto nedostatky odstraˇ nuje Newton˚ uv tvar interpolaˇ cn´ıho polynomu, ve kterém Pn (x) = a0 +a1 (x −x0 )+a2 (x −x0 )(x −x1 )+· · ·+an (x −x0 )(x −x1 ) · · · (x −xn−1 ) . (3.5) Pˇridán´ı dalˇs´ıho uzlu xn+1 je snadné, k Pn (x) staˇc´ı pˇriˇc´ıst jeden dalˇs´ı ˇclen, nebot’ Pn+1 (x) = Pn (x) + an+1 (x − x0 )(x − x1 ) · · · (x − xn ) . Nedostatek (a) Lagrangeova tvaru interpolaˇcn´ıho polynomu jsme tedy pˇrekonali.


Interpolace

Hodnotu z = Pn (¯ x ) urˇc´ıme podobnˇe jako v Hornerovˇe schématu, tj. postupem: z := an a pak pro i = n − 1, n − 2, . . . , 0 poˇc´ıtej z := z(¯ x − xi ) + ai .

(3.6)

Koeficienty ai lze vypoˇc´ıtat pˇr´ımo z interpolaˇcn´ıch podm´ınek (3.1) a dosáhnout tak významné u ´spory v poˇctu operac´ı. Existuje vˇsak jeˇstˇe lepˇs´ı zp˚ usob a ten si ted’ uvedeme. Nejdˇr´ıve definujeme pomˇerné diference : P[xi ] := yi , P[xi , xi+1 ] := (P[xi+1 ] − P[xi ])/(xi+1 − xi ) , P[xi , xi+1 , xi+2 ] := (P[xi+1 , xi+2 ] − P[xi , xi+1 ])/(xi+2 − xi ), a dále, pro 3 ≤ k ≤ n: P[xi , xi+1 , . . . , xi+k−1 , xi+k ] := (P[xi+1 , . . . , xi+k ] − P[xi , . . . , xi+k−1 ])/(xi+k − xi ) . (3.7)


Interpolace

Dá se ukázat, ˇze ai = P[x0 , x1 , . . . , xi ] ,

(3.8)

takˇze Newton˚ uv tvar interpolaˇcn´ıho polynomu je Pn (x) = P[x0 ] + P[x0 , x1 ](x − x0 ) + P[x0 , x1 , x2 ](x − x0 )(x − x1 ) + . . . + P[x0 , x1 , . . . , xn ](x − x0 )(x − x1 ) . . . (x − xn−1 ) . Oznaˇc´ıme-li Pik = P[xi−k , . . . , xi ], pak ai = Pii . Následuje algoritmus v´ ypoˇ ctu pomˇ ern´ ych diferenc´ı Pro i = 0, 1, . . . , n proved’ Pi0 := yi . Pro k = 1, 2, . . . , n provádˇej : pro i = k, k + 1, . . . , n provádˇej : Pik := (Pi,k−1 − Pi−1,k−1 )/(xi − xi−k ) konec cyklu i , konec cyklu k .

(3.9)


Interpolace

x0

P00

x1

P10

P11

x2

P20

P21

P22

x3 .. .

P30 .. .

P31 .. .

P32 .. .

P33

xn

Pn0

Pn1

Pn2

...

..

.

Pn,n−1

Pnn

Výpoˇcet lze pˇrehlednˇe zaznamenat do tabulky, kterou vyplˇ nujeme po sloupc´ıch. Výpoˇcet koeficient˚ u ai = Pii a následný výpoˇcet Pn (¯ x ) podle (3.6) vyˇzaduje 3 1 2 asobic´ıch operac´ı a n2 + 3n operac´ı sˇc´ıtac´ıch. Dosáhli jsme tedy 2 n + 2 n n´ významnˇe niˇzˇs´ıho poˇctu operac´ı, neˇz kolik je jich tˇreba pro výpoˇcet Pn (¯ x) z Lagrangeova interpolaˇcn´ıho polynomu. To znamená, ˇze také nedostatek (b) Lagrangeova interpolaˇcn´ıho polynomu je uspokojivˇe vyˇreˇsen.


Interpolace

Pˇr´ıklad 3.2. Sestav´ıme Newton˚ uv interpolaˇcn´ı polynom pro data z pˇr´ıkladu 3.1. Pr˚ ubˇeh výpoˇctu zaznamenáváme do tabulky. Dostaneme xi

Pi0

Pi1

Pi2

−1

−6

1

−2

2

2

−3

−1

−1

3

2

5

3

Pi3

1

=⇒

a0

=

−6

=⇒

a1

=

2

=⇒

a2

=

−1

=⇒

a3

=

1,

takˇze P3 (x) = −6 + 2 · (x + 1) + (−1) · (x + 1)(x − 1) + 1 · (x + 1)(x − 1)(x − 2) . V bodˇe x¯ = 0,5 podle (3.6) vypoˇcteme P3 (0,5) = ((1 · (0,5 − 2) − 1) · (0,5 − 1) + 2) · (0,5 + 1) − 6 = −1,125 .


Interpolace

Kdyˇz pˇridáme dalˇs´ı uzel x4 = 0 a v nˇem pˇredep´ıˇseme hodnotu y4 = 2, staˇc´ı dopoˇc´ıtat jeden ˇrádek tabulky. Dostaneme x4

P40

P41

P42

P43

P44

0

2

0

2,5

0,5

−0,5

=⇒

a4

a tedy P4 (x) = P3 (x) + (−0,5) · (x + 1)(x − 1)(x − 2)(x − 3) .

=

−0,5 , 2


Interpolace

Chyba aproximace interpolaˇ cn´ım polynomem. Nejdˇr´ıve zavedeme následuj´ıc´ı Oznaˇcen´ı. Symbolem C ha, bi budeme znaˇcit mnoˇzinu vˇsech funkc´ı, které jsou v intervalu ha, bi spojité, a symbolem C k ha, bi pak mnoˇzinu vˇsech funkc´ı, které jsou v intervalu ha, bi spojité spolu se svými derivacemi aˇz do ˇrádu k vˇcetnˇe. Pro k = 0 zˇrejmˇe C 0 ha, bi ≡ C ha, bi. Pˇredpokládejme, ˇze ˇc´ısla yi nejsou libovolná, ale ˇze yi = f (xi ) jsou hodnoty funkce f v uzlech interpolace. Pak nás jistˇe bude zaj´ımat chyba En (¯ x ) := f (¯ x ) − Pn (¯ x) ve zvoleném bodˇe x¯. Pro x¯ = xi je En (xi ) = 0. Jaká je ale chyba mimo uzly interpolace?


Interpolace

Necht’ tedy x¯ je libovolný bod, ha, bi je nˇejaký interval obsahuj´ıc´ı vˇsechny uzly xi interpolace a také zkoumaný bod x¯, a necht’ f ∈ C n+1 ha, bi. Pak pro chybu En (¯ x) plat´ı En (¯ x ) = f (¯ x ) − Pn (¯ x) =

f (n+1) (ξ) (¯ x − x0 )(¯ x − x1 ) . . . (¯ x − xn ) , (n + 1)!

(3.10)

kde ξ = ξ(¯ x ) je (bl´ıˇze neurˇcený) bod z intervalu ha, bi. Zápisem ξ = ξ(¯ x ) pˇritom chceme zd˚ uraznit, ˇze poloha bodu ξ závis´ı nejen na funkci f a na interpolantu Pn , ale také na zvoleném bodu x¯.


Interpolace

Pozn´ amky. Abychom zjednoduˇsili výklad, budeme pˇredpokládat, ˇze x0 < x1 < · · · < xn . 1) Jestliˇze Mn+1 je taková konstanta, ˇze |f (n+1) (x)| ≤ Mn+1 pro kaˇzdé x ∈ ha, bi, pak Mn+1 max |ωn+1 (x)| , |En (¯ x )| ≤ (n + 1)! x∈ha, bi

(3.11)

kde ωn+1 (x) = (x − x0 )(x − x1 ) . . . (x − xn ). Odhad (3.11) je vˇsak obvykle pˇr´ıliˇs pesimistický. 2) Jestliˇze má funkce f (x) derivace vˇsech ˇrád˚ u ohraniˇcené stejnou konstantou, pak pro dostateˇcnˇe velké n je chyba libovolnˇe malá. Pˇr´ıklad 3.3. Pro f (x) = sin x lze vz´ıt Mn+1 = 1, proto |En (x)| ≤

(b − a)n+1 . Dá se dokázat, ˇze (n + 1)!

(b − a)n+1 →0 (n + 1)!

pro

n → ∞,

takˇze Pn (x) → f (x) pro kaˇzdé x z libovolného koneˇcného intervalu ha, bi. 2


Interpolace

1 n=6 n=12 n=18

0.8 0.6 0.4 0.2 0 −1.5

−1

−0.5

0

0.5

Obr. 3.1: Graf funkce |ωn+1 (x)|

1

1.5


Interpolace

3) Jestliˇze interpolaˇcn´ı polynom pouˇz´ıváme k výpoˇctu hodnot interpolované funkce vnˇe intervalu hx0 , xn i, ˇr´ıkáme, ˇze provád´ıme extrapolaci. V tomto pˇr´ıpadˇe m˚ uˇze být chyba aproximace velká, nebot’ hodnota |ωn+1 (x)| rychle roste, kdyˇz se x vzdaluje od x0 doleva nebo od xn doprava. 4) ωn+1 (x) m˚ uˇze nabývat velkých hodnot také uvnitˇr intervalu hx0 , xn i, zejména kdyˇz jsou uzly xi rozm´ıstˇeny rovnomˇernˇe, tj. kdyˇz xi = x0 + ih, kde h je pevnˇe zvolený krok. Na obr. 3.1 vid´ıme, ˇze ve stˇredu intervalu hx0 , xn i nabývá |ωn+1 (x)| nejmenˇs´ıch hodnot, v bl´ızkosti stˇred˚ u krajn´ıch interval˚ u, zejména interval˚ u hx0 , x1 i a hxn−1 , xn i, je vˇsak hodnota |ωn+1 (x)| znaˇcná. Toto chován´ı polynomu ωn+1 (x) se prom´ıtne i do pr˚ ubˇehu interpolantu Pn (x). Pˇr´ıklad 3.4. Sestroj´ıme interpolaˇcn´ı polynom Rungeovy funkce f (x) =

1 1 + 25x 2

na rovnomˇerném dˇelen´ı intervalu h−1, 1i.


Interpolace

2 P10(x) 2

1/(1+25x ) 1.5

1

0.5

0

−0.5 −1

−0.5

0

0.5

1

Obr. 3.2: Interpolaˇcn´ı polynom Rungeovy funkce


Interpolace

Jde o známý pˇr´ıklad, na kterém se demonstruje, ˇze pro rostouc´ı poˇcet d´ılk˚ u chyba interpolace neomezenˇe roste. 2 Vhodným rozm´ıstˇen´ım uzl˚ u lze chybu |ωn+1 (x)| minimalizovat, viz cviˇcen´ı. Tuto moˇznost vˇsak obvykle nemáme, nebot’ uzly jsou pevnˇe dány. Proto pouˇz´ıván´ı interpolaˇcn´ıch polynom˚ u vysokých stupˇ n˚ u obecnˇe nelze doporuˇcit.


Interpolace

Hermitova interpolace. Doposud jsme se zabývali Lagrangeovou interpolac´ı. Jej´ım charakteristickým rysem je to, ˇze interpolaˇcn´ı polynom Pn (x) je urˇcen zadanými hodnotami Pn (xi ) = yi v uzlech xi . Pokud interpolaˇcn´ı polynom urˇcuj´ı nav´ıc také pˇredepsané derivace, hovoˇr´ıme o Hermitovˇe interpolaci. (0) Pˇredpokládejme tedy, ˇze v kaˇzdém uzlu xi je zadáno αi + 1 ˇc´ısel yi , Pn (αi ) (1) yi , . . . , yi . Oznaˇcme α = n + i=0 αi . Pak Hermitovým interpolaˇcn´ım polynomem Pα (x) nazveme polynom stupnˇe nejvýˇse α, který splˇ nuje interpolaˇcn´ı podm´ınky dj (j) Pα (xi ) = yi , dxj

j = 0, 1, . . . , αi ,

i = 0, 1, . . . , n .

(3.12)

Nultou derivac´ı pˇritom rozum´ıme funkˇcn´ı hodnotu. Je dokázáno, ˇze existuje jediný takový polynom. Jestliˇze (j)

yi

=

dj f (xi ) , dxj

j = 0, 1, . . . , αi ,

i = 0, 1, . . . , n ,

ˇr´ıkáme, ˇze Pα (x) je Hermit˚ uv interpolaˇcn´ı polynom funkce f (x).


Interpolace

Necht’ ha, bi je interval obsahuj´ıc´ı uzly interpolace. Jestliˇze f ∈ C α+1 ha, bi, pak pro chybu Hermitovy interpolace v bodˇe x¯ ∈ ha, bi plat´ı f (¯ x ) − Pα (¯ x) =

f (α+1) (ξ) (¯ x − x0 )α0 +1 (¯ x − x1 )α1 +1 . . . (¯ x − xn )αn +1 , (α + 1)!

(3.13)

kde ξ = ξ(¯ x ) je (nˇejaký bl´ıˇze neurˇcený) bod z intervalu ha, bi. Pouˇzit´ı Hermitova polynomu vysokého stupnˇe obecnˇe nelze doporuˇcit, protoˇze (stejnˇe jako u Lagrangeovy interpolace) m˚ uˇze být chyba interpolace mezi uzly znaˇcná. Vzorec pro urˇcen´ı koeficient˚ u Hermitova interpolaˇcn´ıho polynomu je pomˇernˇe komplikovaný, lze ho naj´ıt napˇr. v [Stoer, Bulirsch], zde ho neuvád´ıme. M´ısto toho si na pˇr´ıkladu ukáˇzeme, jak lze Hermit˚ uv polynom urˇcit pˇr´ımo z interpolaˇcn´ıch podm´ınek.


Interpolace

Pˇr´ıklad 3.5. Sestroj´ıme Hermit˚ uv interpolaˇcn´ı polynom pro data podle tabulky xi

yi

yi0

yi00

−1

2

−4

12

1

2

4 (0)

Tedy x0 = −1, α0 = 2, y0 x1 =

(0)

1, α1 = 1, y1

(1)

= 2, y0 (1)

= 2, y1

=

(2)

= −4, y0

= 12,

4.

Protoˇze je pˇredepsáno celkem 5 podm´ınek, Hermit˚ uv polynom navrhneme jako polynom stupnˇe α = 4. Abychom si uˇsetˇrili práci, zap´ıˇseme ho ve tvaru mocninného rozvoje okolo toho bodu, v nˇemˇz je pˇredepsán nejvˇetˇs´ı poˇcet podm´ınek, v naˇsem pˇr´ıpadˇe tedy okolo x0 = −1. Pak P4 (x) = a + b(x + 1) + c(x + 1)2 + d(x + 1)3 + e(x + 1)4 .


Interpolace

Koeficienty a, b, c z´ıskáme snadno. Z podm´ınky P4 (−1) = 2 okamˇzitˇe dostaneme a = 2. Podobnˇe z podm´ınky P40 (−1) = −4 obdrˇz´ıme b = 4 a protoˇze P400 (−1) = 2c, z podm´ınky P400 (−1) = 12 dostaneme c = 6. Dále P4 (1) = 2 − 4 · 2 + 6 · 22 + d · 23 + e · 24 = 2

=⇒

8d + 16e = −16 ,

P40 (1) = −4 + 2 · 6 · 2 + 3 · d · 22 + 4 · e · 23 = 4

=⇒

12d + 32e = −16 .

Tuto soustavu vyˇreˇs´ıme a dostaneme d = −4, e = 1. Tedy P4 (x) = 2 − 4(x + 1) + 6(x + 1)2 − 4(x + 1)3 + (x + 1)4 = x 4 − 1 .

2


Interpolace

Interpolaˇ cn´ı splajny

Jestliˇze chceme interpolovat funkci f (x) na pomˇernˇe dlouhém intervalu ha, bi, mus´ıme ˇzádat splnˇen´ı interpolaˇcn´ıch podm´ınek v dostateˇcnˇe velkém poˇctu bod˚ u. Pokud bude interpolantem polynom, mus´ı být vysokého stupnˇe a to, jak v´ıme, obvykle vede k velkých chybám mezi uzly. Tudy proto cesta nevede. Lepˇs´ı je rozdˇelit interval ha, bi na ˇradu menˇs´ıch subinterval˚ u a na kaˇzdém z nich sestrojit interpolaˇcn´ı polynom niˇzˇs´ıho stupnˇe. Pˇredpokládejme, ˇze a = x0 < x1 < · · · < xi−1 < xi < xi+1 < · · · < xn−1 < xn = b

(3.14)

je dˇelen´ı intervalu ha, bi. V kaˇzdém uzlu xi je pˇredepsána hodnota yi interpolantu. Délku i-tého intervalu hxi−1 , xi i oznaˇc´ıme hi a délku nejdelˇs´ıho intervalu h, tj. hi = xi − xi−1 ,

i = 1, 2, . . . , n ,

h = max hi . 1≤i≤n

(3.15)


Interpolace

Hledaný po ˇcástech polynomický interpolant budeme znaˇcit S(x) a nazveme ho interpolaˇcn´ım splajnem. Na kaˇzdém intervalu hxi−1 , xi i je S(x) polynom, jehoˇz pˇr´ısluˇsnost k i-tému intervalu vyznaˇc´ıme indexem i, tj. S(x) je na intervalu hxi−1 , xi i polynom Si (x) .

K vyjádˇren´ı polynomu Si (x) s výhodou pouˇzijeme lokáln´ı promˇennou s = x − xi−1 . Budeme také pouˇz´ıvat prvn´ı pomˇernou diferenci δi =

yi − yi−1 yi − yi−1 = . xi − xi−1 hi


Interpolace

Line´ arn´ı interpolaˇ cn´ı splajn (dále jen lineárn´ı splajn) je to nejjednoduˇsˇs´ı, co nás napadne: kaˇzdé dva sousedn´ı body [xi−1 , yi−1 ] a [xi , yi ] spoj´ıme u ´seˇckou. Zˇrejmˇe Si (x) = yi−1 +

yi − yi−1 (x − xi−1 ) = yi−1 + sδi xi − xi−1

(3.16)

je lineárn´ı interpolaˇcn´ı polynom procházej´ıc´ı body [xi−1 , yi−1 ] a [xi , yi ]. Lineárn´ı splajn S(x) je spojitá funkce, derivace S 0 (x) je vˇsak ve vnitˇrn´ıch uzlech obecnˇe nespojitá. Jestliˇze yi = f (xi ), i = 0, 1, . . . , n, a f ∈ C 2 ha, bi, pak pro chybu interpolace plat´ı |f (x) − S(x)| ≤ Ch2 , kde x ∈ ha, bi je libovolné a C je konstanta nezávislá na h.

(3.17)


Interpolace

Pro dostateˇcnˇe mnoho uzl˚ u lze uˇcinit chybu libovolnˇe malou. Napˇr´ıklad pˇri vykreslován´ı na obrazovku monitoru s rozliˇsen´ım 1024 x 768 bod˚ u jistˇe staˇc´ı pouˇz´ıt 1024 interpolaˇcn´ıch uzl˚ u k z´ıskán´ı kvalitn´ıho grafu interpolované funkce. Moˇzná bychom byli spokojeni, i kdybychom zvolili ménˇe uzl˚ u, avˇsak pˇri postupném sniˇzován´ı poˇctu uzl˚ u by nutnˇe nastal okamˇzik, kdy by nás jiˇz zaˇcaly ruˇsit ostré hrany grafu S(x) v interpolaˇcn´ıch uzlech. Pokud bychom souˇcasnˇe vykreslovali také funkci f (x), pak by nám zaˇcaly vadit také viditelné odchylky interpolantu S(x) od interpolované funkce f (x) mezi uzly interpolace. Pˇresnˇejˇs´ı interpolant bychom mohli sestrojit tak, ˇze bychom na intervalech hx0 , xk i, hxk , x2k i, . . . aproximovali f (x) pomoc´ı interpolaˇcn´ıch polynom˚ u stupnˇe (nejvýˇse) k > 1. Chyba interpolace by v tom pˇr´ıpadˇe byla u ´mˇerná hk+1 , derivace v uzlech xk , x2k , . . . by vˇsak z˚ ustaly nespojité. Velké k ale nemá smysl pouˇz´ıvat, jinak bychom zase mohli dostat velké chyby mezi uzly interpolace a byli bychom zpˇet v situaci, které jsme se právˇe interpolac´ı po ˇcástech chtˇeli vyhnout.


Interpolace

Velmi populárn´ı je aproximace po ˇcástech kubickým polynomem, která je nejen spojitá, ale má také spojité prvn´ı nebo dokonce i druhé derivace. Popisu takových aproximac´ı se budeme vˇenovat v následuj´ıc´ıch odstavc´ıch. Hermit˚ uv kubick´ y interpolaˇ cn´ı splajn (dále jen Hermit˚ uv kubický splajn) hledáme jako funkci S(x), která a) je v intervalu ha, bi spojitá spolu se svou prvn´ı derivac´ı, tj. S ∈ C 1 ha, bi, b) splˇ nuje interpolaˇcn´ı podm´ınky S(xi ) = yi ,

S 0 (xi ) = di ,

i = 0, 1, . . . , n ,

(3.18)

kde yi , di jsou pˇredepsané funkˇcn´ı hodnoty a derivace, c) je na kaˇzdém intervalu hxi−1 , xi i, i = 1, 2, . . . , n, polynom (nejvýˇse) tˇret´ıho stupnˇe. Si (x) je tedy kubický Hermit˚ uv polynom jednoznaˇcnˇe urˇcený podm´ınkami Si (xi−1 ) = yi−1 , Si0 (xi−1 ) = di−1 , Si (xi ) = yi , Si0 (xi ) = di .


Interpolace

Snadno ovˇeˇr´ıme, ˇze tyto podm´ınky jsou splnˇeny pro Si (x) = yi−1 + sdi−1 + s 2

di−1 − 2δi + di 3δi − 2di−1 − di . + s3 hi hi2

(3.19)

Funkce S(x) je spojitá spolu se svou prvn´ı derivac´ı, druhá derivace uˇz obecnˇe spojitá nen´ı. Jestliˇze yi = f (xi ), di = f 0 (xi ), i = 0, 1, . . . , n, a f ∈ C 4 ha, bi, pak pro chybu interpolace plat´ı |f (x) − S(x)| ≤ Ch4 , (3.20) kde x ∈ ha, bi je libovolné a C je konstanta nezávislá na h. Pokud derivace di nejsou k dispozici, mus´ıme je vypoˇc´ıtat pomoc´ı vhodnˇe zvolených dodateˇcných podm´ınek.


Interpolace

Hermit˚ uv kubick´ y interpolaˇ cn´ı splajn zachov´ avaj´ıc´ı tvar je jednou z moˇznost´ı. Derivace di se vyb´ıraj´ı tak, aby pr˚ ubˇeh S(x) kop´ıroval pr˚ ubˇeh lineárn´ıho splajnu procházej´ıc´ıho body [xi , yi ]. Konkrétnˇe, je-li L(x) lineárn´ı splajn, poˇzadujeme: a) kdyˇz má L(x) ve vnitˇrn´ım uzlu lokáln´ı extrém, pak ho tam má také S(x); b) kdyˇz L(x) mezi dvˇema sousedn´ımi uzly neklesá, pak tam také S(x) neklesá, a podobnˇe, kdyˇz L(x) mezi dvˇema sousedn´ımi uzly neroste, pak tam neroste ani S(x). Jednu zdaˇrilou implementaci lze naj´ıt v MATLABu jako funkci pchip, viz [Matlab], [Moler]. Výpoˇcet smˇernic di prob´ıhá podle následuj´ıc´ıch pravidel: a) Vnitˇrn´ı uzly. Jestliˇze smˇernice δi a δi+1 maj´ı opaˇcná znaménka, nebo je-li nˇekterá z nich rovna nule, tj. kdyˇz plat´ı δi δi+1 ≤ 0 ,

poloˇz´ıme

di = 0 .

V opaˇcném pˇr´ıpadˇe urˇc´ıme di jako zobecnˇený harmonický pr˚ umˇer smˇernic δi a δi+1 ze vztahu w1 w2 w 1 + w2 = + , di δi δi+1

kde w1 = hi + 2hi+1 , w2 = 2hi + hi+1 .

b) Okrajové uzly x0 a xn . Nejjednoduˇsˇs´ı je poloˇzit d0 = δ1 , dn = δn . V algoritmu pchip se ale pouˇz´ıvá kvalitnˇejˇs´ı aproximace vycházej´ıc´ı z kvadratické interpolace, viz [Moler].


Interpolace

Kubick´ y interpolaˇ cn´ı splajn (dále jen kubický splajn). Smˇernice di ve vnitˇrn´ıch uzlech m˚ uˇzeme urˇcit také tak, ˇze poˇzadujeme, aby splajn S(x) mˇel nav´ıc spojitou také druhou derivaci, tj. aby platilo 00 Si00 (xi ) = Si+1 (xi ) ,

i = 1, 2, . . . , n − 1 .

Derivován´ım (3.19) dostaneme Si00 (x) =

(6hi − 12s)δi + (6s − 4hi )di−1 + (6s − 2hi )di . hi2

Pro x = xi je s = hi , takˇze Si00 (xi ) =

−6δi + 2di−1 + 4di . hi

Pro x = xi−1 je s = 0 a Si00 (xi−1 ) =

6δi − 4di−1 − 2di . hi

(3.21)


Interpolace

Kdyˇz v posledn´ım vzorci zvˇetˇs´ıme index i o jedniˇcku, dostaneme 00 Si+1 (xi ) =

6δi+1 − 4di − 2di+1 . hi+1

Dosazen´ım do (3.21) tak dostaneme rovnice hi+1 di−1 + 2(hi+1 + hi )di + hi di+1 = 3(hi+1 δi + hi δi+1 ) ,

i = 1, 2, . . . , n − 1 . (3.22)

Jestliˇze pˇredep´ıˇseme okrajové podm´ınky S 0 (a) = da ,

S 0 (b) = db ,

(3.23)

pak v soustavˇe (3.22) dosad´ıme v prvn´ı rovnici d0 := da a ˇclen h2 da pˇrevedeme na pravou stranu, a v posledn´ı rovnici dosad´ıme dn := db a ˇclen hn−1 db pˇrevedeme na pravou stranu. Soustavu pak ˇreˇs´ıme a z´ıskáme zbývaj´ıc´ı smˇernice di , i = 1, 2, . . . , n − 1.


Interpolace

Matice soustavy je tˇr´ıdiagonáln´ı, diagonálnˇe dominantn´ı, takˇze soustavu lze snadno vyˇreˇsit GEM upravenou pro soustavy s tˇr´ıdiagonáln´ı matic´ı. Jestliˇze yi = f (xi ), i = 0, 1, . . . , n, d0 = f 0 (x0 ), dn = f 0 (xn ), a kdyˇz f ∈ C 4 ha, bi, pak pro chybu interpolace opˇet plat´ı (3.20). Obrázek 3.3 potvrzuje, ˇze pomoc´ı kubického splajnu lze pro data stejná jako v pˇr´ıkladu 3.4 dostat zcela vyhovuj´ıc´ı aproximaci Rungeovy funkce.


Interpolace

Kubický splajn

1

1/(1+25x2) 0.8

0.6

0.4

0.2

0 −1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

Obr. 3.3: Aproximace Rungeovy funkce kubickým splajnem

0.8

1


Interpolace

Kubický splajn má pozoruhodnou extrémáln´ı vlastnost, kterou si ted’ pop´ıˇseme. Oznaˇc´ıme V = {v ∈ C 2 ha, bi| v (xi ) = yi , i = 0, 1, . . . , n, v 0 (x0 ) = d0 , v 0 (xn ) = dn } mnoˇzinu vˇsech funkc´ı, které maj´ı v intervalu ha, bi spojitou druhou derivaci, procházej´ı zadanými body [xi , yi ], i = 0, 1, . . . , n, a v krajn´ıch bodech a = x0 Rb a xn = b jejich derivace nabývaj´ı pˇredepsaných hodnot d0 a dn . Pak a [v 00 (x)]2 dx nabývá na mnoˇzinˇe funkc´ı V své nejmenˇs´ı hodnoty pro kubický splajn S(x), tj. plat´ı Z b Z b 00 2 [S (x)] dx = min [v 00 (x)]2 dx . a

v ∈V

a

Tato vlastnost má zaj´ımavou interpretaci v mechanice. Je totiˇz známo, ˇze ohybová energie homogenn´ıho izotropn´ıho prutu, jehoˇz stˇrednicová ˇcára má Rb rovnici y = v (x), x ∈ ha, bi, má pˇribliˇznˇe hodnotu E (v ) = c a [v 00 (x)]2 dx, kde c je vhodná konstanta. A dále je také známo, ˇze prut, který je donucen procházet pevnými interpolaˇcn´ımi body [xi , yi ] takovým zp˚ usobem, ˇze je namáhán pouze silami p˚ usob´ıc´ımi kolmo k prutu, zaujme pozici s minimáln´ı energi´ı. Extrémáln´ı vlastnost tedy tvrd´ı, ˇze kubický splajn aproximuje stˇrednicovou ˇcáru takového prutu.


Interpolace

Jestliˇze smˇernice da a db v krajn´ıch bodech intervalu ha, bi neznáme, osvˇedˇcil se postup, oznaˇcovaný v anglicky psané literatuˇre jako not a knot. Myˇslenka je jednoduchá: poˇzadujeme, aby splajn byl jednoduchým polynomem tˇret´ıho stupnˇe na prvn´ıch dvou intervalech, tj. pro x0 ≤ x ≤ x2 , a na posledn´ıch dvou intervalech, tj. pro xn−2 ≤ x ≤ xn . V uzlech x1 a xn−1 tedy uˇz nedocház´ı k napojován´ı dvou r˚ uzných polynom˚ u, tj. uzel x1 a xn−1 uˇz nen´ı knot“, ˇcesky uzel splajnu, odtud ” název postupu not a knot“. ” Polynomy S1 (x) a S2 (x) maj´ı v bodˇe x1 spoleˇcnou funkˇcn´ı hodnotu y1 , stejnou prvn´ı derivaci d1 a podle (3.21) také stejnou druhou derivaci. Aby oba polynomy byly totoˇzné staˇc´ı, kdyˇz budou m´ıt v bodˇe x1 také stejnou tˇret´ı derivaci. Stejnou u ´vahu lze provést v bodˇe xn−1 . Dostáváme tak okrajové podm´ınky S1000 (x1 ) = S2000 (x1 ) ,

000 Sn−1 (xn−1 ) = Sn000 (xn−1 ) .

(3.24)


Interpolace

Kdyˇz pomoc´ı (3.19) vyjádˇr´ıme podm´ınku S1000 (x1 ) = S2000 (x1 ) a uprav´ıme ji pomoc´ı prvn´ı rovnice soustavy (3.22), dostaneme rovnici h2 d0 + (h2 + h1 )d1 = [(3h1 + 2h2 )h2 δ1 + h12 δ2 ]/(h1 + h2 ) .

(3.25)

000 Podobnˇe zpracujeme také podm´ınku Sn−1 (xn−1 ) = Sn000 (xn−1 ) a dostaneme rovnici

(hn + hn−1 )dn−1 + hn−1 dn = [hn2 δn−1 + (2hn−1 + 3hn )hn−1 δn ]/(hn−1 + hn ) . (3.26) Neznámé d0 , d1 , . . . , dn pak dostaneme jako ˇreˇsen´ı soustavy rovnic, z nichˇz prvn´ı je rovnice (3.25), pak následuje n − 1 rovnic (3.22) a nakonec pˇrijde jeˇstˇe rovnice (3.26). Pouˇzijeme opˇet GEM upravenou pro soustavy s tˇr´ıdiagonáln´ı matic´ı.


Interpolace

Shrnut´ı. Kubický interpolaˇcn´ı splajn S(x) je funkce, která a) je v intervalu ha, bi spojitá spolu se svou prvn´ı a druhou derivac´ı, tj. S ∈ C 2 ha, bi, b) splˇ nuje interpolaˇcn´ı podm´ınky S(xi ) = yi , i = 0, 1, . . . , n, kde yi jsou pˇredepsané funkˇcn´ı hodnoty, c) je na kaˇzdém intervalu hxi−1 , xi i polynom (nejvýˇse) tˇret´ıho stupnˇe, d) splˇ nuje okrajové podm´ınky (3.23) nebo (3.24).


Interpolace

Pˇr´ıklad 3.6. Oblouk [x(t), y (t)] ≡ [cos t, sin t], t ∈ h0, π/2i, budeme aproximovat kˇrivkou [S x (t), S y (t)], kde S x (t) resp. S y (t) jsou kubické splajny funkc´ı cos t resp. sin t na intervalu h0, π/2i. Interval h0, π/2i rozdˇel´ıme na n stejných d´ılk˚ u, takˇze ti = πi/(2n). Oznaˇc´ıme dix = [S x ]0 (ti ), diy = [S y ]0 (ti ). Protoˇze x 0 (t) = − sin t, y 0 (t) = cos t, poloˇz´ıme d0x = − sin 0 = 0 , dnx = − sin(π/2) = −1 , d0y = cos 0 = 1 , dny = cos(π/2) = 0 . V dalˇs´ım budeme uvaˇzovat n = 3. Soustavu rovnic (3.22) sestav´ıme zvláˇst’ pro x-ovou a zvláˇst’ pro y -ovou sloˇzku. Matice obou soustav je stejná, pravé strany jsou vˇsak r˚ uzné. Na dalˇs´ım ˇrádku je uvedena soustava rovnic se dvˇema pravými stranami a jej´ı ˇreˇsen´ı:   0 1 √ . y 1 1 x  1 4 1 0  d1y = −d2x = 0,865537, −√ 1 2 2 3  d1x d1y  = 3 π . y 1 6 0 1 4 1  d2 d2  d2 = −d1x = 0,499813. − 12 3 2 −1 0 Six (t) a Siy (t) urˇc´ıme podle (3.19). Napˇr´ıklad pro t ∈ hπ/6, π/3i dostaneme . S2x (t) = 0,8660 − 0,4998(t − 61 π) − 0,4431(t − 16 π)2 + 0,1195(t − 61 π)3 , . S2y (t) = 0,5000 + 0,8655(t − 16 π) − 0,2554(t − 16 π)2 − 0,1195(t − 16 π)3 .

2


Interpolace

Konstrukce kubick´ eho interpolaˇ cn´ıho splajnu uˇ zit´ım druh´ ych derivac´ı. Snadno ovˇeˇr´ıme, ˇze kubický polynom Si (x) = yi−1 + s

6δi − 2hi Mi−1 − hi Mi Mi−1 Mi − Mi−1 + s2 + s3 . 6 2 6hi

(3.27)

splˇ nuje podm´ınky Si (xi−1 ) = yi−1 , Si00 (xi−1 ) = Mi−1 , Si (xi ) = yi , Si00 (xi ) = Mi . Funkce S(x), která je na kaˇzdém intervalu hxi−1 , xi i definována pˇredpisem (3.27), proto zˇrejmˇe splˇ nuje podm´ınky S(xi ) = yi , S 00 (xi ) = Mi , i = 0, 1, . . . , n. S(x) je tedy na intervalu ha, bi spojitá a má v nˇem spojitou druhou derivaci. Abychom dostali kubický splajn, mus´ı m´ıt S(x) v ha, bi spojitou také prvn´ı derivaci. Protoˇze nespojitost S 0 (x) m˚ uˇze nastat jedinˇe ve vnitˇrn´ıch uzlech, staˇc´ı poˇzadovat 0 Si0 (xi ) = Si+1 (xi ) ,

i = 1, 2, . . . , n − 1 .

(3.28)


Interpolace

Vyjádˇr´ıme-li (3.28) pomoc´ı (3.27), dostaneme rovnice hi Mi−1 +2(hi +hi+1 )Mi +hi+1 Mi+1 = 6(δi+1 −δi ) ,

i = 1, 2, . . . , n−1 . (3.29)

Kdyˇz zvol´ıme okrajové podm´ınky S 00 (a) = Ma ,

S 00 (b) = Mb ,

(3.30)

dosad´ıme je do (3.29), soustavu rovnic vyˇreˇs´ıme a z´ıskáme Mi , i = 1, 2, . . . , n − 1. Vˇsimnˇete si, ˇze matice soustavy (3.29) je symetrická. Je samozˇrejmˇe moˇzné uvaˇzovat také jiné typy okrajových podm´ınek, napˇr. (3.23) nebo (3.24). Kubický splajn s vlastnost´ı S 00 (a) = S 00 (b) = 0 se nazývá pˇrirozený kubický splajn. Je známo, ˇze pˇrirozený kubický splajn aproximuje pr˚ uhyb prostˇe podepˇreného (homogenn´ıho izotropn´ıho) nosn´ıku namáhaného (silami p˚ usob´ıc´ımi kolmo k nosn´ıku) tak, ˇze tento nosn´ık procház´ı body [xi , yi ]. Mi maj´ı význam ohybových moment˚ u v [xi , yi ].


Interpolace

Interpolace funkc´ı v´ıce promˇ enn´ ych

Omez´ıme se na pˇr´ıpad, kdy f je funkce dvou promˇenných definovaná v oblasti Ω. Interpolace po ˇ c´ astech line´ arn´ı. Pˇredpokládejme, ˇze Ω je mnoho´ uheln´ık. Oblast uheln´ık˚ u T1 , T2 , . . . , Tm , Ω triangulujeme, tj. vyjádˇr´ıme ji jako sjednocen´ı troj´ z nichˇz kaˇzdé dva r˚ uzné bud’to nemaj´ı ˇzádný spoleˇcný bod nebo maj´ı spoleˇcný vrchol popˇr´ıpadˇe maj´ı spoleˇcnou stranu. Mnoˇzinu T = {Tk }m sech takových k=1 vˇ troj´ uheln´ık˚ u nazýváme triangulac´ı oblasti Ω. Vrcholy troj´ uheln´ık˚ u triangulace oznaˇc´ıme P1 = [x1 , y1 ], P2 = [x2 , y2 ], . . . , Pn = [xn , yn ] a nazveme je uzly triangulace. Pˇredpokládejme, ˇze v kaˇzdém uzlu Pi je pˇredepsána hodnota fi = f (xi , yi ) interpolované funkce f .


Interpolace

Po ˇcástech lineárn´ım interpolantem funkce f na oblasti Ω rozum´ıme funkci S, která je v Ω spojitá, splˇ nuje interpolaˇcn´ı podm´ınky S(xi , yi ) = fi , i = 1, 2, . . . , n, a která je na kaˇzdém troj´ uheln´ıku Tk ∈ T lineárn´ı. Na Tk je tedy z = S(x, y ) ≡ Sk (x, y ) rovnice roviny urˇcené hodnotami funkce f ve vrcholech Tk . Pˇripomeˇ nme, ˇze rovnice roviny procházej´ıc´ı body [xa , ya , za ], [xb , yb , zb ] a [xc , yc , zc ] m˚ uˇze být vyjádˇrena ve tvaru x − xa y − ya z − za xb − xa yb − ya zb − za = 0 . xc − xa yc − ya zc − za Vypoˇc´ıtat hodnotu z = S(x, y ) pro (x, y ) ∈ Ω je snadné: urˇc´ıme troj´ uheln´ık Tk , v nˇemˇz bod [x, y ] leˇz´ı, a vypoˇcteme z = Sk (x, y ). Chyba interpolace je t´ım menˇs´ı, ˇc´ım jemnˇejˇs´ı triangulaci zvol´ıme. Kdyˇz f ∈ C 2 (Ω) (tj. kdyˇz f je v Ω spojitá spolu se svými prvn´ımi a druhými parciáln´ımi derivacemi), pak |f (x, y ) − S(x, y )| ≤ Ch2 , kde h je nejdelˇs´ı strana troj´ uheln´ık˚ u triangulace a C je konstanta nezávislá na h.


Interpolace

Interpolace po ˇ c´ astech biline´ arn´ı. Pˇredpokládejme, ˇze Ω = ha, bi × hc, di je obdéln´ık. Pomoc´ı dˇelen´ı a = x0 < x1 < · · · < xn = b, c = y0 < y1 < · · · < ym = d rozloˇz´ıme výchoz´ı obdéln´ık Ω na menˇs´ı obdéln´ıky Rij = {(x, y ) | xi−1 ≤ x ≤ xi , yj−1 ≤ y ≤ yj }, i = 1, 2, . . . , n, j = 1, 2, . . . , m. Pˇredpokládejme, ˇze v uzlech [xi , yj ] jsou pˇredepsány hodnoty fij = f (xi , yj ) funkce f. Na obdéln´ıku Rij definujeme funkci Sij (x, y ) = fi−1,j−1

yj − y x − xi−1 yj − y xi − x + fi,j−1 + xi − xi−1 yj − yj−1 xi − xi−1 yj − yj−1

+ fi−1,j

xi − x y − yj−1 x − xi−1 y − yj−1 + fij . xi − xi−1 yj − yj−1 xi − xi−1 yj − yj−1

Funkce Sij je bilineárn´ı , tj. pro pevné x = C je Sij (C , y ) lineárn´ı funkce promˇenné y a pro pevné y = D je Sij (x, D) lineárn´ı funkce promˇenné x. Funkce Sij je interpolant funkce f na obdéln´ıku Rij , tj. Sij nabývá ve vrcholech [xi−1 , yj−1 ], [xi , yj−1 ], [xi−1 , yj ] a [xi , yj ] stejných hodnot jako funkce f , jak se snadno pˇresvˇedˇc´ıme.


Interpolace

Na Ω definujeme funkci S pˇredpisem S(x, y ) = Sij (x, y ) pro (x, y ) ∈ Rij , i = 1, 2, . . . , n, j = 1, 2, . . . , m. Funkce S je po ˇcástech ( = po obdéln´ıc´ıch Rij ) bilineárn´ı. Protoˇze S(xi , yj ) = fij , i = 0, 1, . . . , n, j = 0, 1, . . . , m, ˇrekneme, ˇze S je po ˇcástech bilineárn´ı interpolant funkce f na obdéln´ıku Ω. Snadno ovˇeˇr´ıme, ˇze S je v Ω spojitá (staˇc´ı si uvˇedomit, ˇze Sij , a tedy také S, je na kaˇzdé stranˇe obdéln´ıka Rij jednoznaˇcnˇe urˇcena pomoc´ı hodnot funkce f v koncových bodech této strany). Kdyˇz f ∈ C 2 (Ω), pak pro chybu interpolace plat´ı odhad |S(x, y )−f (x, y )| ≤ C (h2 +k 2 ) , a C je konstanta nezávislá na h, k.

kde h = max (xi −xi−1 ) , k = max (yj −yj−1 ) 1≤i≤n

1≤j≤m


Obsah

3



Metoda nejmenˇs´ıch ˇ ctverc˚ u



oznaˇcuje postup pro pˇribliˇzné ˇreˇsen´ı pˇreurˇcených nebo nepˇresnˇe zadaných soustav rovnic, zaloˇzený na minimalizaci kvadrát˚ u jejich rezidu´ı. Prokl´ ad´ an´ı dat kˇrivkami je významná skupina u ´loh, které lze metodou nejmenˇs´ıch ˇctverc˚ u ˇreˇsit (v anglicky psané literatuˇre se pro tyto aplikace pouˇz´ıvá ´lohách jde. oznaˇcen´ı curve fitting). Popiˇsme si, o co v takových u Necht’ tedy t je nezávisle promˇenná, napˇr´ıklad ˇcas, a y (t) je neznámá funkce promˇenné t, kterou chceme aproximovat. Pˇredpokládejme, ˇze jsme provedli m pozorován´ı, tj. hodnoty y byly zmˇeˇreny pro urˇcité (navzájem r˚ uzné) hodnoty t: yi = y (ti ) ,

i = 1, 2, . . . , m .

Naˇs´ım zámˇerem je modelovat y (t) lineárn´ı kombinac´ı n bázových funkc´ı pro nˇejaké n ≤ m: y (t) ≈ x1 ϕ1 (t) + x2 ϕ2 (t) + · · · + xn ϕn (t) := Rn (t) . Bázové funkce navrhujeme podle oˇcekávaného pr˚ ubˇehu neznámé funkce y (t), urˇcit se maj´ı parametry x1 , x2 , . . . , xn . Funkce Rn (t) se ve statistice nazývá lineárn´ı regresn´ı funkce.



Návrhová matice A modelu je obdéln´ıková matice, která má m ˇrádk˚ u a n sloupc˚ u:   ϕ1 (t1 ) ϕ2 (t1 ) . . . ϕn (t1 )    ϕ1 (t2 ) ϕ2 (t2 ) . . . ϕn (t2 )     A=  .. .. ..  ≡ (ϕ1 , ϕ2 , . . . , ϕn ) ,  . . .    ϕ1 (tm ) ϕ2 (tm ) . . . ϕn (tm ) kde ϕi = (ϕi (t1 ), ϕi (t2 ), . . . ϕi (tm ))T je i-tý sloupec A. Maticová formulace modelu je y ≈ Ax , kde y = (y1 , y2 , . . . , ym )T jsou namˇeˇrená data a x = (x1 , x2 , . . . , xn )T je vektor neznámých parametr˚ u. Symbol ≈ vyjadˇruje pˇribliˇznou rovnost.



y

i 2 i

r

y

i+1

2 i+1

r =y −R (t ) i

i

n

r

i

2

ri−1 y

i−1

r2m

ym 2

r2

y2 y1

Σm r2 → min i=1 i

r21 t1

t2

ti−1

ti

Obr. 3.4: Princip metody nejmenˇs´ıch ˇctverc˚ u

ti+1

tm



Rezidua jsou rozd´ıly mezi pozorován´ımi a modelem: ri = yi − Rn (ti ) = yi −

n X

ϕj (ti )xj ≡ yi −

j=1

n X

aij xj ,

i = 1, 2, . . . , n ,

j=1

kde aij = ϕj (ti ). V maticovém zápisu r = y − Ax .

(3.31)

Parametry xi chceme urˇcit tak, aby rezidua byla co nejmenˇs´ı. Metodu nejmenˇs´ıch ˇctverc˚ u dostaneme, kdyˇz minimalizujeme souˇcet ˇctverc˚ u rezidu´ı: krk2 :=

m X i=1

ri2 → min .

(3.32)



Nˇekdy se pouˇz´ıvá také váˇzená metoda nejmenˇs´ıch ˇctverc˚ u: kdyˇz jsou nˇekterá pozorován´ı d˚ uleˇzitˇejˇs´ı nebo pˇresnˇejˇs´ı neˇz ostatn´ı, pak m˚ uˇzeme pˇrisoudit jednotlivým pozorován´ım r˚ uzné váhy wi > 0 a minimalizovat souˇcet váˇzených ˇctverc˚ u m X 2 krkw := wi ri2 → min . i=1

Je-li napˇr´ıklad chyba i-tého pozorován´ı pˇribliˇznˇe rovna ei , zvol´ıme wi = 1/ei . Je dobré si uvˇedomit, ˇze kaˇzdou metodu pro ˇreˇsen´ı neváˇzené metody nejmenˇs´ıch ˇctverc˚ u lze pouˇz´ıt pro ˇreˇsen´ı metody váˇzené: staˇc´ı pronásobit yi a i-tý ˇrádek A √ souˇcinitelem wi .



Pˇribliˇzné ˇreˇsen´ı pˇreurˇcené soustavy rovnic Ax = y (tj. rovnic je v´ıc neˇz neznámých), které minimalizuje délku rezidua r = y − Ax, se nazývá ˇreˇsen´ı soustavy lineárn´ıch rovnic metodou nejmenˇs´ıch ˇctverc˚ u. ˇ sen´ı minimalizaˇcn´ı u ´lohy (3.32) mus´ı splˇ novat Norm´ aln´ı soustava rovnic. Reˇ nutnou podm´ınku pro extrém: 2  m n 2 X X ∂krk ∂ yi − k = 1, 2, . . . , n . = aij xj  = 0 , ∂xk ∂xk i=1

j=1

Kdyˇz provedeme naznaˇcené derivován´ı, dostaneme   m n X X ∂krk2 yi − =2 aij xj  (−aik ) = 0 ∂xk i=1

a odtud

n m X X j=1

! aik aij

i=1

xj =

j=1

m X

aik yi ,

k = 1, 2, . . . , n ,

i=1

coˇz lze zapsat maticovˇe jako AT Ax = AT y .

(3.33)



Soustava lineárn´ıch rovnic (3.33) je známa jako normáln´ı soustava rovnic. Kdyˇz jsou sloupce matice A lineárnˇe nezávislé, je matice G := AT A pozitivnˇe definitn´ı ´lohy (3.32), tj. plat´ı a ˇreˇsen´ı x∗ normáln´ı soustavy rovnic je jediné ˇreˇsen´ı u ky − Ax∗ k2 = min ky − Axk2 . x

Vyjádˇr´ıme-li normáln´ı soustavu rovnic pomoc´ı vektor˚ u ϕi , dostaneme (ϕ1 , ϕ1 )  (ϕ , ϕ )  2 1  ..   . (ϕn , ϕ1 ) 

    (ϕ1 , ϕ2 ) . . . (ϕ1 , ϕn ) x1 (ϕ1 , y)     (ϕ2 , ϕ2 ) . . . (ϕ2 , ϕn )    x2   (ϕ2 , y)    =   .. ..   ..   ..  ,  .   .  . ... . (ϕn , ϕ2 ) (ϕn , ϕn ) xn (ϕn , y)

kde (ϕk , ϕj ) =

m X i=1

ϕk (ti )ϕj (ti )

a

(ϕk , y) =

m X

(3.34)

ϕk (ti )yi

i=1

jsou skalárn´ı souˇciny vektor˚ u ϕk , ϕj a ϕk , y. Matice G soustavy (3.34) se nazývá Gramova matice soustavy vektor˚ u ϕj , j = 1, 2, . . . , n.



Pˇri návrhu aproximace Rn (t) bychom mˇeli vyb´ırat funkce ϕi (t) tak, aby sloupce ϕi matice A byly lineárnˇe nezávislé. V opaˇcném pˇr´ıpadˇe, jak se dá ukázat, má u ´loha (3.32) nekoneˇcnˇe mnoho ˇreˇsen´ı, coˇz je zˇrejmˇe neˇzádouc´ı. Uved’me si dva významné speciáln´ı pˇr´ıpady, pro které jsou sloupce matice A lineárnˇe nezávislé (d˚ ukaz lze naj´ıt napˇr. v [Berezin]): a) ϕj (t) je polynom stupnˇe j − 1, napˇr. ϕj (t) = t j−1 , j = 1, 2, . . . , n; b) pro n = 2N + 1, kde N je celé nezáporné ˇc´ıslo, vol´ıme ϕ1 (t) = 1,

ϕ2k (t) = cos kt,

ϕ2k+1 (t) = sin kt,

k = 1, 2, . . . , N ,

a ˇcasy pozorován´ı“ ti vyb´ıráme z intervaluhc, c + 2π), kde c je libovolné ” ˇc´ıslo. Aproximace Rn (t) je v pˇr´ıpadˇe a) algebraický polynom a v pˇr´ıpadˇe b) trigonometrický polynom.



Kdyˇz je m = n a matice A je regulárn´ı, pak x∗ = A−1 y a r = o, tj. Rn (ti ) = yi , i = 1, 2, . . . , m. Pokud jsou vˇsak namˇeˇrená data yi zat´ıˇzena chybami, pak nen´ı u ´ˇcelné, aby funkce Rn (t) tyto chyby kop´ırovala. Naopak, aby Rn (t) vˇerohodnˇe vystihovala (rekonstruovala) neznámou funkci y (t), je ˇzádouc´ı, aby Rn (t) namˇeˇrená data vyrovnávala (vyhlazovala). To je ale moˇzné jen tehdy, kdyˇz poˇcet pozorován´ı m je výraznˇe vˇetˇs´ı neˇz poˇcet n návrhových parametr˚ u, tj. pro m n.



Pˇr´ıklad 3.7. Pro data pˇredepsaná tabulkou ti

0

0,5

1

1,5

2

2,5

3

yi

3,57

2,99

2,62

2,33

2,22

2,10

2,05

urˇc´ıme aproximaci R2 (t) = x1 + x2 e −t metodou nejmenˇs´ıch ˇctverc˚ u. Zˇrejmˇe ϕ1 (t) = 1 a ϕ2 (t) = e −t . Normáln´ı soustava rovnic je tvaru     7  7 7 P P P 1 · yi x1 1 · e −ti 1·1      i=1  i=1 i=1 .   =        P 7 7 7 P −ti P −ti −ti −ti e · yi x2 e ·1 e ·e i=1

i=1

i=1

Vypoˇcteme-li pˇr´ısluˇsné sumy, dostaneme soustavu (zobrazujeme nejvýˇse 4 desetinná m´ısta) ! ! ! 7 2,4647 x1 17,88 = , 2,4647 1,5805 x2 7,4422 . . jej´ıˇz ˇreˇsen´ı je x1 = 1,9879 a x2 = 1,6087. Hledaná aproximace . . −t R2 (t) = 1,99 + 1,61e a krk = 0,0651. 2



Pˇr´ıklad 3.8. Daty z pˇredchoz´ıho pˇr´ıkladu proloˇz´ıme postupnˇe polynomy prvn´ıho, druhého a tˇret´ıho stupnˇe. Za bázové vol´ıme funkce ϕj (t) = t j−1 , kde j = 1, 2, . . . , n a n = 2, 3, 4. a) Polynom prvn´ıho stupnˇe. Pro ϕ1 (t) = 1 a ϕ2 (t) = t dostaneme normáln´ı soustavu rovnic ! ! ! 7 10,5 x1 17,88 = , 10,5 22,75 x2 23,45 . . . která má ˇreˇsen´ı x1 = 3,28 a x2 = −0,48, takˇze R2 (t) = 3,28 − 0,48t . a krk = 0,4756. Aproximace lineárn´ım polynomem tedy nen´ı vhodná, nebot’ je málo pˇresná. b) Polynom druhého stupnˇe. Normáln´ı soustava rovnic      7 10,5 22,75 x1 17,88       10,5     22,75 55,125     x2  =  23,45  22,75 55,125 142,1875 x3 49,0625 . . . má ˇreˇsen´ı x1 = 3,53, x2 = −1,09 a x3 = 0,20, takˇze . . R3 (t) = 3,53 − 1,09t + 0,2t 2 a krk = 0,1006. Velikost rezidua se zmenˇsila, je vˇsak stále vˇetˇs´ı neˇz v pˇr´ıkladu 3.7.



c) Polynom tˇret´ıho stupnˇe. Normáln´ı soustava rovnic (zobrazujeme nejvýˇse 4 desetinná m´ısta)      7 10,5 22,75 55,125 x1 17,88       10,5     22,75 55,125 142,1875     x2   23,45     =    22,75   x3  49,0625 55,125 142,1875 381,2813      55,125

142,1875

381,2813

1049,5469

x4

116,78

. . . . má ˇreˇsen´ı x1 = 3,57, x2 = −1,35, x3 = 0,43 a x4 = −0,05, takˇze . . R4 (t) = 3,57 − 1,35t + 0,43t 2 − 0,05t 3 a krk = 0,0360.



Pokud bychom stupeˇ n polynomu dále zvyˇsovali, zjistili bychom, ˇze polynom R7 (t) ˇsestého stupnˇe procház´ı vˇsemi body [ti , yi ], takˇze je to interpolaˇcn´ı polynom. Vˇsimnˇete si jeˇstˇe prvk˚ u Gramových matic: s rostouc´ım ˇrádem nejvˇetˇs´ı koeficient prudce roste. Spolu s n´ım prudce roste také ˇc´ıslo podm´ınˇenosti Gramových matic. Do následuj´ıc´ı tabulky jsme zaznamenali ˇc´ısla podm´ınˇenosti κ2 (G) pro n = 2, 3, . . . , 7 (spoˇctená v MATLABu pomoc´ı maticové k · k2 normy) n

2

3

4

5

6

7

κ2 (G)

16

4,27 · 102

1,91 · 104

1,20 · 106

1,17 · 108

2,31 · 1010

Pokud bychom tabulku dat P z pˇr´ıkladu 3.7 rozˇs´ıˇrili o dalˇs´ı sloupce, mohli bychom n teoreticky v regresn´ı funkci j=1 xj t j−1 zvˇetˇsovat n (aˇz do poˇctu m sloupc˚ u). Prakticky to vˇsak moˇzné nen´ı, nebot’ pro velké n by ˇc´ıslo podm´ınˇenosti Gramovy matice bylo ne´ unosnˇe velké. Regresn´ı funkce Rn (t) s bázovými funkcemi ϕj (t) = t j−1 , j = 1, 2, . . . , n, se proto pro vˇetˇs´ı n nepouˇz´ıvaj´ı. Pokud potˇrebujeme polynomickou regresn´ı funkci vyˇsˇs´ıho stupnˇe, mˇeli bychom pouˇz´ıt tzv. ortogonáln´ı polynomy, viz napˇr. [Stoer, Bulirsch]. 2



ˇ ıslo podm´ınˇenosti κ2 (A) je definováno ˇ sen´ı pˇreurˇ Reˇ cen´ ych soustav rovnic. C´ také pro obdéln´ıkovou matici A a plat´ı κ2 (AT A) = [κ2 (A)]2 . To znamená, ˇze ˇc´ıslo podm´ınˇenosti Gramovy matice AT A je výraznˇe vˇetˇs´ı neˇz ˇc´ıslo podm´ınˇenosti návrhové matice A. Normáln´ı rovnice se proto k ˇreˇsen´ı rozsáhlých soustav pˇreurˇcených rovnic nehod´ı. Existuj´ı lepˇs´ı postupy vyuˇz´ıvaj´ıc´ı tzv. QR rozklad matice A nebo tzv. pseudoinverzi matice A,viz napˇr. [Stoer, Bulirsch], [Moler]. Obˇe tyto metody jsou stabiln´ı a Gramovu matici nepouˇz´ıvaj´ı. Kvalitn´ı implementace jsou dostupné napˇr. v MATLABu, viz [Matlab]. Tyto metody se také bez pot´ıˇz´ı vypoˇrádaj´ı s pˇr´ıpadem, kdyˇz jsou sloupce matice A lineárnˇe závislé. Dá se ukázat, ˇze v tom pˇr´ıpadˇe má pˇreurˇcená soustava rovnic nekoneˇcnˇe mnoho ˇreˇsen´ı (ve smyslu metody nejmenˇs´ıch ˇctverc˚ u, tj. má nekoneˇcnˇe mnoho ˇreˇsen´ı, která maj´ı stejnou délku rezidua). Jestliˇze A+ je pseudoinverze matice A, pak A+ y je jednoznaˇcnˇe definované ˇreˇsen´ı, které minimalizuje délku rezidua, a pokud je takových ˇreˇsen´ı nekoneˇcnˇe mnoho, pak A+ y je to z nich, které má nejmenˇs´ı délku.


Obsah

3



Literatura


Literatura

ˇ ˇ I.S. Berezin, N.P. Zidkov: Cislennyje metody I,II, Nauka, Moskva, 1962. G. Dahlquist, G. ˚ A Bj˝ork: Numerical Methods, Prentice-Hall, Englewood Cliffs, NJ, 1974. M. Fiedler: Speciáln´ı matice a jejich pouˇzit´ı v numerické matematice, SNTL, Praha, 1981. D. Hanselman, B. Littlefield: Mastering MATLAB 7, Pearson Prentice Hall, Upper Saddle River, NJ, 2005. G. H¨ ammerlin, K. H. Hoffmann: Numerical Mathematics, Springer-Verlag, Berlin, 1991. M. T. Heath: Scientific Computing. An Introductory Survey, McGraw-Hill, New York, 2002. I. Horov´ a, J. Zelinka: Numerické metody, uˇcebn´ı text Masarykovy univerzity, Brno, 2004. J. Kobza: Splajny, uˇcebn´ı text Palackého univerzity, Olomouc, 1993. J. Klapka, J. Dvoˇrák, P. Popela: Metody operaˇcn´ıho výzkumu, uˇcebn´ı text, FSI VUT Brno, 2001.


Literatura

J. H. Mathews, K. D. Fink: Numerical Methods Using MATLAB, Pearson Prentice Hall, New Jersey, 2004. MATLAB: Mathematics, Version 7, The MathWorks, Inc., Natick, 2004. G. Meurant: Computer Solution of Large Linear Systems, Elsevier, Amsterodam, 1999. S. M´ıka: Numerické metody algebry, SNTL, Praha, 1985.. C. B. Moler: Numerical Computing with MATLAB, Siam, Philadelphia, 2004. http://www.mathworks.com/moler. J. Nocedal, S. J. Wright: Numerical Optimization, Springer Series in Operations Research, Springer, Berlin, 1999. A. Quarteroni, R. Sacco, F. Saleri: Numerical Mathematics, Springer, Berlin, 2000. W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling: Numerical Recipes in Pascal, The Art of Scientific Computing, Cambridge University Press, Cambridge, 1990.


Literatura

P. Pˇrikryl: Numerické metody matematické analýzy, SNTL, Praha, 1985. A. R. Ralston: Základy numerické matematiky, Academia, Praha, 1973. K. Rektorys: Pˇrehled uˇzité matematiky I,II, Prometheus, Praha, 1995. J. Stoer, R. Bulirsch: Introduction to Numerical Analysis, Springer-Verlag, New York, 1993. E. Vitásek: Numerické metody, SNTL, Praha, 1987. W. Y. Yang, W. Cao, T. S. Chung, J. Morris: Applied Numerical Methods Using Matlab, John Willey & Sons, New Jersey, 2005.

Aproximace funkcí. Numerické metody. Ústav matematiky. 1. listopadu 2006

Recommend Documents