Interpolace a aproximace Interpolace algebraickým polynomem a p g ý p y aproximace metodou nejmenších čtverců
Interpolace algebraickým polynomem
Aproximace metodou nejmenších čtverců
Úloha. Dána tabulka hodnot x i, yi, xi≠xj pro i ≠j. Hodnoty jsou „přesné“. Hledáme polynom nejvýše n polynom nejvýše n‐tého tého stupně, p(x)=a0+a1x+a2x2+…+anxn (tj. koeficienty a0.. an), který prochází všemi zadanými body tj p(x i)=y i. všemi zadanými body, tj. p(x
Úloha. Dána tabulka hodnot xi, yi, kde x i může být = xj pro i ≠j. Hodnoty jsou zatíženy chybou, např. experimentálně naměřené . Hledáme takový polynom nejvýše n‐tého ýp y j ý stupně, p(x)=a0+a1x+a2x2+…+anxn (tj. koeficienty a0.. an), aby součet druhých mocnin odchylek p(x oc odc y e p( i)) od y od yi by byl minimální. á
Rovnice pro určení koeficientů hl d éh hledaného polynomu sestavujeme z l j podmínky p(x i)=y i. Stupeň hledaného polynomu je určen počtem zadaných bodů. ů
Rovnice pro určení koeficientů hledaného polynomu sestavujeme z podmínky minimality polynomu sestavujeme z podmínky minimality kvadratické odchylky . Stupeň hledaného polynomu je určen předpokládanou závislostí hodnot předpokládanou závislostí hodnot.
Interpolace algebraickým polynomem. Úloha. Dána tabulka hodnot xi, yi, xi≠xj pro i ≠j. x
x0
x1
x2
…
xn
y
y0
y1
y2
…
yn
Hledáme polynom nejvýše n‐tého stupně, p(x)=a0+a1x+a2x2+…+anxn (tj. koeficienty a0.. an), který prochází všemi zadanými body, tj. p(xi)=yi. Rovnice pro určení koeficientů hledaného polynomu sestavujeme z podmínky p(xi)=yi. ⎛1 x0 x02 L x0n ⎞⎛ a0 ⎞ ⎛ y0 ⎞ bod [x 0, yy 0]: a ]: a0+a1x 0 +a2x 0 2+…+a + +anx 0 n = y = y0 ⎜ ⎟⎜ ⎟ ⎜ ⎟ 2 n 2 n bod [x 1, y 1]: a0+a1x 1 +a2x 1 +…+anx 1 = y 1 ⎜1 x1 x1 L x1 ⎟⎜ a1 ⎟ ⎜ yn ⎟ = ⎜ ⎟ bod [x 2, y 2]: a0+a1x 2 +a2x 2 2+…+anx 2 n = y 2 ⎜M M ⎟ ⎜ ⎟ M M ⎟ M ⎜ ⎜ ⎟ ⎜ ⎟ … 2 n ⎟⎜ ⎟ ⎜y ⎟ ⎜ a L 1 x x x 2 n n⎠ n n n ⎠⎝ n ⎠ bod [x n, y n]: a0+a1x n +a2x n +…+anx n = y n ⎝144 123 ⎝{ 42 4443 zadané x -ové hodnoty
neznámé
Jsou li všechny x‐ové hodnoty navzájem různé, je matice soustavy regulární, tedy existuje jediné řešení (a0,a1,a2,…,+an)T , tj. právě jeden polynom p(x). Stupeň hledaného polynomu je určen počtem bodů. Stupeň hledaného polynomu je určen počtem bodů
zadané
Interpolace algebraickým polynomem ‐ příklad Dána tabulka hodnot :
x
‐2
‐1
0
1
y
30
11
4
9
a) Sestavte soustavu lineárních rovnic pro určení interpolačního polynomu. b) Určete interpolační polynom. c) Vypočítejte interpolovanou hodnotu v bodě x=0.5. d) Výsledky zobrazte. a)) Hledáme polynom polynom p y p y nejvýše 3. stupně: p(x)=a j ý p p( ) 0+a1x+a2x2+a3x3 tj. neznámé j koeficienty a0, a1, a2, a3 . Soustavu rovnic získáme dosazením zadaného bodu do p(x). bod [‐2 , 30]: p(‐2) = a0+a1(‐2)+a2(‐2)2+a3(‐2)3= 30 ⎫ ⎛1 − 2 4 − 8 ⎞⎛ a0 ⎞ ⎛ 30 ⎞ ⎜ ⎟⎜ ⎟ ⎜ ⎟ bod [ 1 11]: p( 1) = a0+a1(‐1)+a bod [‐1 , 11]: p(‐1) = a ( 1)+a2(‐1) ( 1)2+a3(‐1) ( 1)3= 11 ⎪ ⎪ ⎜1 − 1 1 − 1 ⎟⎜ a1 ⎟ ⎜ 11 ⎟ ⇒ = bod [0 , 4]: p(0) = a0+a1(0)+a2(0)2+a3(0)3 = 4 ⎬ ⎜1 0 0 0 ⎟⎜ a ⎟ ⎜ 4 ⎟ ⎪ ⎜ ⎟⎜ 2 ⎟ ⎜ ⎟ 2 3 bod [1 , 9]: p(1) = a0+a1(1)+a2(1) +a3(1) = 9 ⎪ ⎜1 1 1 1 ⎟⎜ a ⎟ ⎜ 9 ⎟
⎭
⎝
b) Vyřešíme soustavu: a0=4, a1=‐1, a2 =6, a3 =0, ttím jsme určili polynom 2. stupně: p(x)=4 – js e u č po y o . stup ě: p( ) x + 6x 6 2 c) Interpolovaná hodnota v bodě x = 0.5 je p(0.5) = 4 ‐ 0.5 + 6(0.25) = 5
⎠⎝
3
⎠ ⎝
⎠
Interpolace algebraickým polynomem ‐ příklad d) Zobrazení výsledků:
Interpolace polynomem ‐ ověření v MATLABu ověření v MATLABu xii = [‐2 ‐1 0 1]; [2 0 ] % d é i % zadané xi yi = [30 11 4 9]; %zadané yi x=0.5; % zadaný bod x n=size(xi); % n(1) : počet řádků, n(2): počet sloupců xi A=[ones(n)’ cumprod(xi’*ones(1,n(2)‐1),2] % matice soustavy %%nebo “ručně“ ručně , neefektivně: neefektivně: A=[ones(n 1) xi’ (xi A=[ones(n,1) xi (xi’)).^2 (xi ^2 (xi’)).^3] ^3] ai = A\yi % řešení soustavy, koeficienty ai p = fliplr(ai’) % „otočení“ ai, y = polyval(p, x) % výpočet hodnoty pro x = 0.5 %%%%%%%%%%% zobrazení výsledků %%%%%%%%% xx = linspace(xi(1), xi(n(2)), 20) linspace(xi(1), xi(n(2)), 20) % rozdělíme interval <‐2,1> % rozdělíme interval 2,1 na 20 bodů na 20 bodů yy = polyval(p, xx) % vypočteme yy = p(xx) plot(xx,yy,xi,yi, ’ks’, x, y, ’g*’) % zobraz
Aproximace metodou nejmenších čtverců Úloha. Dána tabulka hodnot xi, yi, kde x i může být = xj pro i ≠j. Hledáme takový polynom nejvýše n‐tého stupně, p(x)=a Hledáme takový polynom nejvýše n‐tého stupně p(x)=a0+a1x+a2x2+…+a + +anxn (tj. koeficienty a0… an), aby součet druhých mocnin odchylek p(xi) od yi byl minimální.
Rovnice pro určení koeficientů hledaného polynomu sestavujeme z podmínky n
δ = ∑ ( p( xi ) − yi ) → min minimalityy druhé mocniny kvadratické odchylky : D = . y y y 2
2
i =0
Minimum kvadratické odchylky je v bodech, kde parciální derivace podle neznámých a0… aan jjsou nulové.
Nalezený polynom má ze všech polynomů stejného stupně nejmenší kvadratickou Nalezený polynom má ze všech polynomů stejného stupně nejmenší kvadratickou odchylku, tedy nejlépe ve smyslu metody nejmenších čtverců aproximuje zadaná data.
Odvození soustavy rovnic pro polynom 2. stupně p(x)=a0+a1x+a2x2 ,x=x0…xn,y=y0…yn: n
Z íš Zapíšeme druhou mocninu kvadratické odchylky: D= d h i k d ti ké d h lk D δ = ∑ ( p ( xi ) − yi ) 2 2
i =0
D = (p(x0)‐y0)2 + (p(x1)‐y1)2 + (p(x2)‐y2)2 …+(p(xn)‐yn)2 tj. D = (a j ( 0+a1x0+a2x02‐yy0) 2+ (a ( 0+a1x1+a2x12‐yy1) 2 +(a ( 0+a1x2+a2x22‐yy2) 2 … +(a ( 0+a1xn+a2xn2‐yyn) 2 Vyjádříme parciální derivace D a položíme je = 0. ∂D = 2(a0 + a1 x0 + a2 x02 − y0 ) + 2(a0 + a1 x1 + a2 x12 − y1 ) + 2(a0 + a1 x2 + a2 x22 − y2 ) + L + 2(a0 + a1 xn + a2 xn2 − yn ) = 0 ∂a0
∂D = 2(a0 + a1 x0 + a2 x02 − y0 ) x0 + 2(a0 + a1 x1 + a2 x12 − y1 ) x1 + 2(a0 + a1 x2 + a2 x22 − y2 ) x2 + L + 2(a0 + a1 xn + a2 xn2 − yn ) xn = 0 ∂a1 ∂D = 2(a0 + a1 x0 + a2 x02 − y0 ) x02 + 2(a0 + a1 x1 + a2 x12 − y1 ) x12 + 2(a0 + a1 x2 + a2 x22 − y2 ) x22 + L + 2(a0 + a1 xn + a2 xn2 − yn ) xn2 = 0 ∂a2
Rovnice upravíme a soustavu rovnic zapíšeme maticově. (a0 + a0 + a0 + L + a0 ) + (a1 x0 + a1 x1 + a1 x2 + L + a1 xn ) + (a2 x02 + a2 x12 + a2 x22 + L + a2 xn2 ) = y0 + y1 + y2 + L + yn (a0 x0 + a0 x1 + a0 x2 + L + a0 xn ) + (a1 x02 + a1 x12 + a1 x22 + L + a1 xn2 ) + (a2 x03 + a2 x13 + a2 x23 + L + a2 xn3 ) = x0 y0 + x1 y1 + L + xn yn (a0 x02 + a0 x12 + a0 x22 + L + a0 xn2 ) + (a1 x03 + a1 x13 + a1 x23 + L + a1 xn3 ) + (a2 x04 + a2 x14 + a2 x24 + L + a2 xn4 ) = x02 y0 + x12 y1 + x22 y2 + L + xn2 yn ⎫ a0 (n + 1) + a1 ∑ xi + a2 ∑ x = ∑ yi ⎪ ⎛⎜ n + 1 i =0 i =0 i =0 ⎪ ⎜ n n n n ⎪ ⎜ n a0 ∑ xi + a1 ∑ xi2 + a2 ∑ xi2 = ∑ xi yi ⎬ ⇒ ⎜ ∑ xi i =0 i =0 i =0 i =0 ⎪ ⎜ i =n0 n n n n ⎪ ⎜ 2 a0 ∑ xi2 + a1 ∑ xi3 + a2 ∑ xi4 = ∑ xi2 yi ⎪ ⎜ ∑ xi i =0 i =0 i =0 i =0 ⎭ ⎝ i =0 n
n
2 i
n
n
∑x i =0 n
i
∑x i =0 n
2 i
∑x i =0
3 i
⎞ ⎞ ⎛ n x ⎟ ⎜ ∑ yi ⎟ ∑ i =0 ⎟⎛ a0 ⎞ ⎜ i =0 ⎟ n n ⎜ ⎟ ⎟ ⎜ ⎟ xi3 ⎟⎜ a1 ⎟ = ⎜ ∑ xi yi ⎟ ∑ i =0 ⎟⎜ a2 ⎟ ⎜ i =n0 ⎟ n ⎝ ⎠ ⎜ 4⎟ 2 ⎟ xi ⎟ ∑ ⎜ ∑ xi yi ⎟ i =0 ⎠ ⎝ i =0 ⎠ n
2 i
Aproximace metodou nejmenších čtverců ‐ příklad. Dána tabulka hodnot:
x
‐2
‐1
‐1
0
0
4
y
‐7
‐1.8
‐2.2
1.1
0.9
‐7
a) Určete polynom nejvýše 2. stupně, který ve smyslu metody nejmenších čtverců nejlépe aproximuje data daná tabulkou hodnot. b) Vypočítejte kvadratickou odchylku. Vypočítejte kvadratickou odchylku c) Výsledky zobrazte. a) Hledáme p(x)=a0+a1x+a2x2 . Neznámé a0, a1, a2 určíme ze soustavy rovnic: kde : n + 1 = 6 (počet bodů v tabulce) n n ⎛ ⎛ n ⎞ 2⎞ ⎜ n +1 ⎜ ⎜ n ⎜ ∑ xi ⎜ i =n0 ⎜ x2 ⎜∑ i ⎝ i =0
∑x ∑x i =0 n
i
∑x i =0 n
2 i
∑ xi3 i =0
⎟ ⎜ ∑ yi ⎟ i =0 ⎟⎛ a0 ⎞ ⎜ i =0 ⎟ n n ⎜ ⎟ ⎟ ⎜ ⎟ xi3 ⎟⎜ a1 ⎟ = ⎜ ∑ xi yi ⎟ ∑ i =0 ⎟ ⎟⎜ a2 ⎟ ⎜ i =n0 n ⎝ ⎠ ⎜ 4⎟ 2 ⎟ xi ⎟ ∑ ⎜ ∑ xi yi ⎟ i =0 ⎠ ⎝ i =0 ⎠ i
5
5
∑x i =0
i
5
∑y i =0
i
5
5
= 0, ∑ x = 22, ∑ x = 54, ∑ xi4 = 274 i =0
2 i
i =0
3 i
i =0
5
5
i =0
i =0
= −16, ∑ xi yi = −10, ∑ xi2 yi = −144
2. ⎛ 6 0 22 ⎞⎛ a0 ⎞ ⎛ − 16 ⎞ ⎛ a0 ⎞ ⎛ 1 ⎞ Tedy: p(x)=1+2x ‐ x ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ 0 22 54 = − 10 a ⎜ ⎟⎜ 1 ⎟ ⎜ ⎟ ⇒ ⎜ a1 ⎟ = ⎜ 2 ⎟ ⎜ 22 54 274 ⎟⎜ a ⎟ ⎜ − 144 ⎟ ⎜ a ⎟ ⎜ − 1⎟ ⎝ ⎠⎝ 2 ⎠ ⎝ ⎠ ⎝ 2⎠ ⎝ ⎠
Aproximace metodou nejmenších čtverců ‐ příklad.
b) K zadané tabulce přidáme hodnoty nalezeného polynomu v zadaných bodech, K zadané tabulce přidáme hodnoty nalezeného polynomu v zadaných bodech vypočteme odchylky p(x)‐y a jejich druhé mocniny. Součet čísel v posledním řádku je druhá mocnina kvadratické odchylky : δ2. Výslednou δ získáme odmocněním. x
‐2 2
‐1 1
‐1 1
0
0
4
y
‐7
‐1.8
‐2.2
1.1
0.9
‐7
p( ) p(x)
‐7 7
‐2 2
‐2 2
1
1
‐7 7
(p(x)‐y)
0
‐0.2
0.2
‐0.1
0.1
0
(p(x)‐y) (p(x) y)2
0
0 04 0.04
0 04 0.04
0 01 0.01
0 01 0.01
0
c) Zobrazení
p(x)=1+2x ‐ x2 δ2 = 0.1, δ = 0.3162
Aproximace metodou nejmenších čtverců – v MATLABu. xi = [‐2 ‐1 ‐1 0 0 4]; yi = [‐7 ‐1.8 ‐2.2 1.1 0.9]; N i ( i) N=size(xi); A=[N(2) sum(xi) sum(xi.^2); sum(xi) sum(xi.^2) sum(xi.^3); sum(xi.^2) sum(xi.^3) sum(xi.^4)] B=[sum(yi); sum(xi.*yi); sum((xi.^2).*yi)] ai = A\B; = A\B; poly = fliplr(ai) ; P = polyval(poly, xi) D = P – yi; Delta2=sum(D.^2); delta = sqrt(Delta2); %%%% zobraz xx=linspace(xi(1), xi(N(2)),20); yy=polyval(poly, xx); p ( , yy, , y , plot(xx, yy, xi, yi, ’k*’); );
