Mat-1.414, V2 kevät 2000
Luennot viikolla 5
- 1
- Interpolaatio 2
- ...
- n
x=linspace(-1,1);n=10;
kertomat=gamma(1:n+1);
for j=1:n+1
P(:,j)=polyval(fliplr(1./kertomat(1:j)),x)';
end;
plot(x,P)
hold on
plot(x,exp(x))
plot(x,exp(x),x,P(:,3))
>> figure
>> plot(x,exp(x),x,P(:,3))
>> plot(x,exp(x),x,P(:,3),x,P(:,5))
>> plot(x,exp(x)'-P(:,3))
>> plot(x,exp(x)'-P(:,5))
>> plot(x,exp(x)'-P(:,6)) % Tarkkaile mittakaavan muuttumista
>> plot(x,exp(x)'-P(:,7))
>> plot(x,exp(x)'-P(:,8))
>> plot(x,exp(x)'-P(:,9))
>> plot(x,exp(x)'-P(:,10))
>>
Taulukointi on tosi kätevää
Kannattaa yleensä katsoa vain osaa, esim. 7 ensimmäistä riviä ja kaikki
sarakkeet:
>> P(1:7,:)
ans =
Columns 1 through 7
1.0000 0 0.5000 0.3333 0.3750 0.3667 0.3681
1.0000 0.0202 0.5002 0.3434 0.3818 0.3743 0.3755
1.0000 0.0404 0.5008 0.3535 0.3889 0.3821 0.3832
1.0000 0.0606 0.5018 0.3637 0.3961 0.3900 0.3910
1.0000 0.0808 0.5033 0.3738 0.4036 0.3981 0.3989
1.0000 0.1010 0.5051 0.3840 0.4112 0.4063 0.4071
1.0000 0.1212 0.5073 0.3942 0.4191 0.4147 0.4154
Columns 8 through 11
0.3679 0.3679 0.3679 0.3679
0.3754 0.3754 0.3754 0.3754
0.3830 0.3830 0.3830 0.3830
0.3909 0.3909 0.3909 0.3909
0.3988 0.3988 0.3988 0.3988
0.4070 0.4070 0.4070 0.4070
0.4153 0.4153 0.4153 0.4153
Kuinka monta termiä otettava, jotta |virhe| < tol
Lasketaan e:lle likiarvo kaavalla
> e=Sum(1/k!,k=0..n);
n
-----
\ 1
e = ) ----
/ k!
-----
k = 0
tol=10^(-6)
N=0:14
gamma(N+1)> 3/tol % Muista: gamma(n+1)=n!
>> N(ans)
ans =
10 11 12 13 14
Siis 10 termiä riittää.
Annettu dataa: (esim. n=7)
> matrix([[seq(x[i],i=0..7)],[seq(y[i],i=0..7)]]);
[x[0] x[1] x[2] x[3] x[4] x[5] x[6] x[7]]
[ ]
[y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7]]
Lagrangen interpolaatiopolynomit, Lagrangen kertojafunktiot ("Cardinal functions") etc. ...
Tässä eräs Matlab-toteutus:
function y=lagrange(xdata,xev,j)
% Lagrangen kertojafunktio ("Cardinal function") HA maalisk. -95
% y=lagrange(xdata,xeval,j);
% Lasketaan xdata:aan liittyv{ j:s l-funktio pisteiss{ xeval
% esim: y=lagrange([1 2 3 4],[1 2 3 4],0) antaa 1 0 0 0
%
j=j+1; % Matlab-vektorissa indeksi alkaa 1:st{
y=1; tuloind=1:length(xdata);ind=tuloind~=j;tuloind=tuloind(ind);
for i=tuloind;y=y.*(xev-xdata(i))./(xdata(j)-xdata(i));end;
Interpolaatiopolynomin P(x) arvo lasketaan näin:
n
-----
\
P(x) = ) y[i] L[i](x)
/
-----
i = 0
Matlabissa halutaan tietysti laskea x-vektorissa. No, mikä ettei:
Tehdään matriisi:
Lmat= [Lagrange(xd,x,0)',[Lagrange(xd,x,1)',...,[Lagrange(xd,x,n)'];
y=L*yd' % Kätevää!!
Maplella katsottu ln-esimerkki
Kun datapisteitä on vain 3, tehdään käsin, ilman eleganttia matriisikertolaskua.
xdata=[9.0,9.5,11.0];ydata=log(xdata);
x=linspace(8,12);
L0=lagrange(xdata,x,0);
L1=lagrange(xdata,x,1);
L2=lagrange(xdata,x,2);
p=ydata(1)*L0+ydata(2)*L1+ydata(3)*L2;
plot(x,p,x,log(x),xdata,ydata,'o')
plot(x,p-log(x),xdata,zeros(size(xdata)),'o')
No toimiikohan tuo matriisikertolasku tosiaan? Katsotaan!
y=[L0' L1' L2']*ydata';
plot(x,y,x,log(x),xdata,ydata,'o')
Jehuu!
Interpolaatio valmiilla funktioilla, polyfit/polyval
Tästähän jo oli täällä
Luentotehtävä
Tee edellinen ln-esimerkki polyfit/polyval-tyylillä.