Polynomanpassning i MATLAB

(1)

ACPU - 29 January 2004 Yngve Sundblad Yngve Sundblad Föreläsning 4 sid.1 SF 1518/19 ht 2015 9 sept.

Polynomanpassning i MATLAB

Funktionsanropet c=polyfit(x,y,n) ger koefficiemterna i ett n:e-gradspolynom som anpassar sig till y-värdena för  

x-värdena med lämplig metod. 

I tredje föreläsningens exempel 3, där andragradspolynomet  y = F(x) = c₁ + c₂x + c₃x² anpassas till 

x=[15:1:19]';y=[15 8 5 3 3]'  ger c=polyfit(x,y,2) 

[1.0714 -39.3286 363.6000]’ 

I exemplet A\y ger [363.6 -39.3286 1.0714]

yv=polyval(c,xv) ger värdet av polynomet för xv.

polyval(c,17) ger 4.6571,  

stämmer med 4.6571-2.9000(x-17)+1.0714(x-17)^2 (OBS! för lab 1:5)

OBS! Omvänd ordning, högsta koefficienten först, också i c för  polyval 

Frågor från F3: norm(x,Inf) och norm(x,inf) är samma sak  Finns ingen funktion för att sortera om till diagonaltung (god idé)

(2)

Interpolation med polynom

Interpolation innebär att beräkna en polynomkurva som går genom ett antal punkter och använda detta för att uppskatta mellanliggande värden. (Man kan använda andra basfunktioner men det ligger utanför kursen.)

Extrapolation innebär att använda polynomkurvan för att uppskatta värden utanför de givna.

Algebrans fundamentalsats (Gauss 1799):

Givet n punkter (värdepar) (x₁, y₁), (x₂, y₁),…(x_n, y_n). 

Det finns precis ett polynom P av grad n-1 som uppfyller P(x_k)=y_k

 ”Naiva” formen: P(x) = c₁+ c₂x + c₃x²+ … + c_nx^n-1, 

kräver någon form av ekvationslösning för att beräkna c:na Newtons interpolationsformel (1675): 

P(x) = c₁+ c₂(x-x₁) + c₃(x-x₁)(x-x₂) + … + c_n(x-x₁)(x-x₂)…(x-x_n-1)

(3)

Beräkning av koefficienterna

Med ”naiva” blir systemmatrisen och ekvationen 

( )( )=( )

A*c = y, i MATLAB c = A\y Möjliga problem:

Stora x_k ger stora element i A, blir illa konditionerat, dvs små fel i y_k blir stora fel i c_k. Lösning: Centrera kring x_medel

Stort gradtal, n, ger stora svängningar (oberoende av metod)

1 x₁x₁² … x₁^n-1 1 x₂x₂² … x₂^n-1  

…

1 x_nx_n² … x_n^n-1

c₁ c₂  

… c_n

y₁ y₂  

… y_n

(4)

Interpolationsproblemen i MATLAB

c=polyfit(x,y,n-1) kan också användas för att beräkna interpolationspolynomets genom n punkter koefficienter, liksom yv=polyval(c,xv) som ger värdet av polynomet för xv 

OBS! Ordningen, högsta koefficienten först.

Exempel: Följande data anger andelen internetanvändare i

procent i hela världen olika år, anpassa med sjättegradspolynom T=[1997 2000 2003 2006 2009 2012 2015]; 

a = [2 7 12 18 26 36 42]; 

c6=polyfit(T,a,6);

Warning: Polynomial is badly conditioned. Add points with distinct X values,

reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.  

> In polyfit (line 75) c6 =1.0e+15 *

0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0062 2.0606

(5)

Centrera, först, minska gradtal sedan

T=T-2006; format long; c6=polyfit(T,a,6) 

-0.000007620789514 -0.000102880658436 -0.000171467764060 0.00771604938271 0.113271604938272 2.272222222222221 17.999999999999996 (I short nollor i c(1))

format short; 

a6_2018=polyval(c6,12); ger 23.0000

% Pröva med varannan punkt och tredjegradare 

T=[-9:6:9]; a=[2 12 26 42]; c3=polyfit(T,a,3) 

-0.0015 0.0417 2.3472 18.6250 

a3_2018=polyval(c3,12); ger 50.1250 Tpl=-9:0.1:9;  

in6pl=[c6(1)*Tpl.^6+c6(2)*Tpl.^5+c6(3)*Tpl.^4+c6(4)*T pl.^3+c6(5)*Tpl.^2+c6(6)*Tpl+c6(7)];

in4pl=[c4(1)*Tpl.^3+c4(2)*Tpl.^2+c4(3)*Tpl+c4(4)];

plot(2006+T,a,’*’); hold on; 

plot(2006+Tpl,in6pl,2006+Tpl,in3pl);

(6)

Interpolationspolynomen

Den röda.

6:egrads, verkar böja, rentav neråt. 

Värdet 23   år 2018   bekräftar. 

Den blå, 3:egrads,   verkar   snällare.

Värdet 50

bekräftar.

(7)

Runges fenomen 1

Ännu värre problem kan inträffa i utkanterna av en interpolation. I

Pohl,sid.139, interpoleras e^x/(1+9x²) i tretton punkter i intervallet (-1,1) med polynom av tolfte graden:

x=-1:1/6:1; y=exp(x)./(1+9*x.^2); 

c=polyfit(x,y,12)

170.5044 -19.6793 -449.9423 51.9314 445.0745 -51.3695 -211.9794 24.4662 53.6316 -6.1891 -8.1345 0.9578 1.0000

xpl=-1:0.01:1;

ypl=zeros(1,201);

for i=1:13

ypl=ypl+c(i)*xpl.^(13-i);

end

plot(xpl, ypl, xpl, exp(xpl)./(9+xpl.^2))

(8)

Runges (Carl, 1856-1927) fenomen 2

Den röda är  funktionen,  den blå   polynomet  som  

interpolerar 

där det ska  

men ändå  

svänger  

våldsamt 

i kanterna!

(9)

Styckvis interpolation

Man kan få ”snällare” interpolationskurvor genom att interpolera styckvis mellan varje par av punkter och byta polynom i varje punkt. Extra bra blir det om derivatan (lutningen) och till och med andraderivatan (krökningen) blir oförändrad vid punkterna.

Enklast är styckvis linjär interpolation, dvs med linjer mellan punkterna som följer på varandra.

P(t) = P(x_k+t*h_k)=y_k+t*∆y_k, där h_k=x_k+1-x_k och ∆y_k=y_k+1-y_k, 0≤t≤1 I allmänhet blir det hörn (diskontinuerlig derivata) i punkterna, olämpligt t.ex. för profiler som ska ha lågt luftmostånd.

Om punkterna ligger tätt kan man lura ögat, som vid plottning i MATLAB.

Styckvis kvadratisk används sällan, just inte bättre, men styckvis kubisk interpolation används mycket.

(10)

Hermiteinterpolation (1864)

Mellan varje par av x lägger man ett tredjegradspolynom, vars värden resp. derivator är desamma i ändpunkterna som de angänsandes. Därmed försvinner hörnen.

Tredjegradspolynomet P_k(x) mellan x_k och x_k+1 ska alltså ha egenskaperna

P_i(x_k)=y_k, P_k(x_k+1)=y_k+1,P_k’(x_k)=y’(x_k)=d_k, P_k’(x_k+1)=y’(x_k+1)=d_k+1 Vi ansätter ett tredjegradspolynom på hanterbar form:

P_k(x)=c₁ + c₂(x-x_k)+c₃(x-x_k)²(x-x_k+1)+c₄(x-x_k)(x-x_k+1)²

Genom att derivera, sätta in x_k och x_k+1,använda h_k=x_k+1-x_k och

∆y_k=y_k+1-y_kfår vi

c₁=y_k, c₂=∆y_k/h_k,c₃=(d_k+1-c₂)/h_k², c₄=(d_k-c₂)/h_k² Beteckningar som i Pohl, dock .._k istf .._id istf k

(11)

Hermite, exempel

Exempel (Pohl 4.13):

x = [1 5 9]; y= [2.54 6.52 6.02]; yprim=[2.12 0.22 -0.33];

Två polynom, P₁ i intervallet 1≤x≤5, P₂ i intervallet 5≤x≤9 P₁: c₁=y₁=2.54, c₂=∆y₁/h_1,=(6.52-2.54)/4=0.995

c₃=(d₂-c₂)/h₁² =(0.22-0.995)/16, c₄=(d₁-c₂)/h₁²=(2.12-0.995)/16 P₁(x) =2.54+0.995(x-1)+(x-1)(x-5)(1.125(x-5)-0.775(x-1))/16 Använd MATLABs symboliska algebra:

syms x;  

P1=2.54+0.995*(x-1)+

(x-1)*(x-5)*(1.125*(x-5)-0.775*(x-1))/16;

ger (199*x)/200 + (((7*x)/20 - 97/20)*(x - 1)*(x - 5))/16 + 309/200

P1=simplify(P1)

ger (7*x^3)/320 - (139*x^2)/320 + (4677*x)/1600 + 47/1600

(12)

Hermite, exempel, forts

P₂: c₁=y₂=6.52, c₂=∆y₂/h_2,=(6.02-6.52)/4=-0.125

c₃=(d₃-c₂)/h₂² =(-0.33+0.125)/16, c₄=(d₂-c₂)/h₂²=(0.22+0.125)/16 P₂(x) =6.52-0.125(x-5)+(x-5)(x-9)(0.345(x-9)-0.205(x-5))/16

Använd MATLABs symboliska algebra:

P2=simplify(6.52-0.125*(x-5)+(x-5)*(x-9)* 

(0.345*(x-9)-0.205*(x-5))/16) ger (7*x^3)/800 - (101*x^2)/400  

+ (1671*x)/800 + 259/200

ezplot(P1,[1,9]) ger den blå  hold on 

ezplot(P2,[1,9]) ger den röda Kurvorna går mjukt över i varandra   för x=5 (mitten)

(13)

Kubiska splines

I skepps-, flyg-, tåg- och bilbranschen behöver man mjuka kurvor där också andraderivatorna är kontinuerliga.

Detta löste man med elastiska kurvlinjaler (eng. spline, svenska ri) som med god noggrannhet visar sig styckvis likna

tredjegradspolynom.

Med Hermite-interpolation med derivatorna   ersatta: y’(x_k)= d_k = ∆y_k/h_k)=d_k och  

andraderivatorna kontinuerliga och   andraderivatan=0 i kanterna får man  ett tridiagonalt ekvationssystem, som  ger värdena till en spline (se Pohl 4.14) 

(14)

Hermite och splines i MATLAB

I MATLAB finns funktionen

pchip - Piecewise Cubic Hermite Interpolating Polynomial  som ger yi-värdena för x-värdena i xi 

om den anropas för givna x- och y-vektorer med   Herm = pchip(x,y,X), 

där Herm sätts till interpolerade värdena för X-vektorn och för funktionen spline 

Ri = spline(x, y, X) 

ger Ri värdena för X-vektorn av den splinekurva   som definieras av x och y

(15)

Överkurs

I NAM (Eriksson), avsnitt 3.4-3.7 beskrivs för interpolation 3.4 Hackigheter i kurvan genom illa-konditionering

3.5 Experimentell störningsanalys

3.6 Polynommultiplikation genom interpolation 3.7 Flerdimensionell interpolation

Det ingå inte i kursen och kommer inte på tentan, men är läsvärt för speciellt intresserade.

Polynomanpassning i MATLAB