888888888888 88 88
88 88 ,d ,d88
88 88 88 888888
88 ,adPPYba, 88,dPPYba, MM88MMM 88
88 a8P_____88 88P' "8a 88 88
88 8PP""""""" 88 88 88 88
88 "8b, ,aa 88 88 88, 888 88
88 `"Ybbd8"' 88 88 "Y888 888 88
%bernstein
x=linspace(5,17); %halutaan polynomeja välille [5,17]
a=5;
b=17;
p=(x-a)/(b-a); %kuvataan väli [5,17] välille [0,1]
b1=bernstein(5,4,0); %muodostetaan bernsteinin polynomi välille [0,1]
y1=polyval(b1,p); %lasketaan sen arvoja
plot(x,y1) %piirretään arvot alkuperäisen välin [5,17] muuttujan x funktiona
[a,b]=ginput; %ammutaan kuvasta sopiva kohta tekstille: klikkaus hiirellä ja enter
text(a,b,'B^5_4') %laitetaan ammuttuun kohtaan teksti
b2=bernstein(7,3,0); %toistetaan
y2=polyval(b2,p);
hold on
plot(x,y2)
[a,b]=ginput;
text(a,b,'B^7_3')
b3=bernstein(9,2,0); %toistetaan
y3=polyval(b3,p);
plot(x,y3,'k')
[a,b]=ginput;
text(a,b,'B^9_2')
print -djpeg 'h6t1bernsteineja.jpg' %tallennetaan kuva
%legendre
le1=legepo(7,0); %käytetään funktiota legepo, tehdään ensimmäinen polynomi
z1=polyval(le1,p); %tehdään samat piirtotyöt kuin aiemmin
plot(x,z1)
[a,b]=ginput;
text(a,b,'L_7')
le2=legepo(3,0); %toistetaan
z2=polyval(le2,p);
hold on
plot(x,z2,'k')
[a,b]=ginput;
text(a,b,'L_3')
print -djpeg 'h6t1legendreja.jpg' %tallennetaan kuva
le3=lege01(4,0); %käytetään kokeeksi toista funktiota lege01
z3=polyval(le3,p);
plot(x,z3,'k')
[a,b]=ginput;
text(a,b,'L_4')
hold on
le4=lege01(5,0); %toistetaan
z4=polyval(le4,p);
plot(x,z4,'m')
[a,b]=ginput;
text(a,b,'L_5')
print -djpeg 'h6t1legendreja2.jpg' %tallennetaan kuva
888888888888 88 ad888888b,
88 88 ,d d8" "88
88 88 88 a8P
88 ,adPPYba, 88,dPPYba, MM88MMM ,d8P"
88 a8P_____88 88P' "8a 88 a8P"
88 8PP""""""" 88 88 88 a8P'
88 "8b, ,aa 88 88 88, 888 d8"
88 `"Ybbd8"' 88 88 "Y888 888 88888888888
%funktio f on tallennettu nimellä h6t2
%harjoituksissa 5 käytiin läpi koodipätkä
n=2;
x=linspace(0,1);
F=zeros(1,n+1);
for k=0:n
F=F+h5t3(k/n)*bernstein(n,k,n+1);
end
y=polyval(F,x);
f=h5t3(x);
plot(x,f) %piirretään f sinisellä
hold on
plot(x,y,'r') %piirretään F punaisella
for k=0:n
K(k+1,:)=[k/n,h5t3(k/n)];
end
K=[K;K(1,:)];
plot(K(:,1),K(:,2),'k')
legend('funktio','approksimaatio','rajoittava kolmio')
title('Harjoitus 5, tehtävä 3')
%poimitaan siitä tarpeelliset osat ja muutetaan jotain
%tehdään aluksi b2 ja b32 ja kuvat niistä f:n kanssa
n=2;
x=linspace(0,1);
F=zeros(1,n+1);
for k=0:n
F=F+h6t2(k/n)*bernstein(n,k,n+1);
end
eval(['b' num2str(n) '=F;']) %poimitaan polynomi tässä vaiheessa nimelle bn, missä n on luku
y=polyval(F,x); %piirretään malliksi kuva
f=h6t2(x);
plot(x,f) %piirretään f sinisellä
hold on
plot(x,y,'r') %piirretään b2 punaisella
eval(['legend(''f'',''b' num2str(n) ''')'])
title('Harjoitus 6, tehtävä 2')
print -djpeg 'h6t2_f_ja_b2.jpg' %tallennetaan kuva
n=32;
x=linspace(0,1);
F=zeros(1,n+1);
for k=0:n
F=F+h6t2(k/n)*bernstein(n,k,n+1);
end
eval(['b' num2str(n) '=F;']) %poimitaan polynomi tässä vaiheessa nimelle bn, missä n on luku
y=polyval(F,x); %piirretään malliksi kuva
f=h6t2(x);
plot(x,f) %piirretään f sinisellä
hold on
plot(x,y,'r') %piirretään b2 punaisella
eval(['legend(''f'',''b' num2str(n) ''')'])
title('Harjoitus 6, tehtävä 2')
print -djpeg 'h6t2_f_ja_b32.jpg' %tallennetaan kuva
%tehdään sitten polynomeja silmukalla
for m=1:10
n=2^m;
x=linspace(0,1);
F=zeros(1,n+1);
for k=0:n
F=F+h6t2(k/n)*bernstein(n,k,n+1);
end
eval(['b' num2str(n) '=F;']) %poimitaan polynomi tässä vaiheessa nimelle bn, missä n on luku
end
%nyt on saatu polynomit b2,b4,...,b1024
%jos oltaisiin laitettu m=1:20, niin olisi tullut liian pitkä silmukka ja b2^20 olisi ollut liian paha laskettavaksi
%tällöin silmukka olisi saatu rikki komennolla Ctrl+c
%lasketaan virheitä saaduille polynomeille
virhe(b2,100) %lasketaan polynomin b2 virhe käyttämällä 100:aa evaluointipistettä
ans =
1.7744
virhe(b2,10000)
ans =
1.7745
%10000:lla pisteellä saadaan suurinpiirtein sama arvo virheelle
%katsotaan virheet b2,...,b1024
eva=10000;
for m=1:10
e(m)=eval(['virhe(b' num2str(2^m) ',' num2str(eva) ');']);
end
e'
ans =
1.0e+043 *
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
4.703241552379511
NaN
NaN
NaN
%tuli hurjia lukuja
%katsotaan viisi ensimmäistä
e(1:5)'
ans =
1.774520397313124
1.088464445617766
0.772270037293300
0.494781504366073
0.293193400609904
%siis virhe pienenee
%tutkitaan komentoa loglog
eva=10^3;
for m=1:10
n=10*m;
F=zeros(1,n+1);
for k=0:n
F=F+h6t2(k/n)*bernstein(n,k,n+1);
end
e(m)=virhe(F,eva);N(m)=n;
end
loglog(N,e)
print -djpeg 'h6t2virhe.jpg'
%nyt kuvassa näkyy ehkä suora
888888888888 88 ad888888b,
88 88 ,d d8" "88
88 88 88 a8P
88 ,adPPYba, 88,dPPYba, MM88MMM aad8"
88 a8P_____88 88P' "8a 88 ""Y8,
88 8PP""""""" 88 88 88 "8b
88 "8b, ,aa 88 88 88, 888 Y8, a88
88 `"Ybbd8"' 88 88 "Y888 888 "Y888888P'
%tehdään f:n approksimaatio legendren polynomeilla kannassa P_5
n=5;
le=zeros(1,n+1);
for k=0:n
le=le+lege01(k,n+1)*kerroin(k);
end
x=linspace(0,1);
y=h6t2(x);
z=polyval(le,x);
plot(x,y)
hold on
plot(x,z,'k')
eval(['legend(''f'',''approks. kannassa P_' num2str(n) ''')'])
title('Harjoitus 6, tehtävä 3')
%katsotaan saadaanko Matlabista virheilmoitus
kerroin(1000)
Warning: Infinite or Not-a-Number function value encountered.
> In quadl at 107
In kerroin at 4
ans =
NaN
%saatiin virheilmoitus
%kuinka iso virhe on?
virhe(le,100)
ans =
0.4634
%tämä on hieman pienempi kuin polynomin b16 virhe. siis nyt 5. asteen legendren polynomi
%voittaa 16. asteen bernsteinin polynomin maksimivirhessä
%tehdään approksimaatiot le1,...,le20 ja piirretään kuva niiden virheistä
for m=1:20
n=m;
le=zeros(1,n+1);
for k=0:n
le=le+lege01(k,n+1)*kerroin(k);
end
e(m)=virhe(le,10^4);
end
plot(1:20,e)
xlabel('n')
ylabel('approksimaation le_n virhe')
title('Harjoitus 6, tehtävä 3, virheen pieneneminen')
print -djpeg 'h6t3_virhe.jpg'
plot(1:20,log(e))
xlabel('n')
ylabel('approksimaation le_n virheen logaritmi')
title('Harjoitus 6, tehtävä 3, virheen pieneneminen2')
print -djpeg 'h6t3_virhe2.jpg'
%kuvan h6t3_virhe2.jpg mukaan virhe pienenee
888888888888 88 ,d8
88 88 ,d ,d888
88 88 88 ,d8" 88
88 ,adPPYba, 88,dPPYba, MM88MMM ,d8" 88
88 a8P_____88 88P' "8a 88 ,d8" 88
88 8PP""""""" 88 88 88 8888888888888
88 "8b, ,aa 88 88 88, 888 88
88 `"Ybbd8"' 88 88 "Y888 888 88
%jos otetaan n=k-1, niin interpolointipolynomi kulkee datapisteiden kautta
k=10; %valitaan datapisteiden lukumäärä
n=k-1; %valitaan iterpolaatiopolynomin aste
d=linspace(0,1,k); %valitaan datapisteet tasavälisesti
D=h6t2(d); %lasketaan funktion f arvo niissä, saadaan tason pisteet (d,D=f(d))
eval(['la' num2str(n) '=polyfit(d,D,n);']) %muodostetaan Lagrangen interpolointipolynomi la_n
x=linspace(0,1); %valmistaudutaan piirtämään
y=h6t2(x); %lasketaan funktion arvoja piirtopisteissä
eval(['z=polyval(la' num2str(n) ',x);']) %lasketaan polynomin la_n arvo niissä
plot(x,y) %piirretään kuvat
hold on %piirretään lisää
plot(x,z,'k')
eval(['legend(''f'',''Lagrangen interp. polynomi la' num2str(n) ''')'])
title('Harjoitus 6, tehtävä 4')
%miltä virheet näyttävät?
virhe(la9,10^3)
ans =
0.1394
>> virhe(la4,10^3)
ans =
0.7751
%tehdään silmukalla erilaisia interpolointipolynomeja ja katsotaan virhettä
for m=1:10
k=m; %valitaan datapisteiden lukumäärä
n=k-1; %valitaan iterpolaatiopolynomin aste
d=linspace(0,1,k); %valitaan datapisteet tasavälisesti
D=h6t2(d); %lasketaan funktion f arvo niissä, saadaan tason pisteet (d,D=f(d))
la=polyfit(d,D,n); %muodostetaan Lagrangen interpolointipolynomi la_n
e(m)=virhe(la,10^3);
end
plot(1:10,e)
title('Harjoitus 6, tehtävä 4, virhe')
xlabel('n')
ylabel('approksimaation la_n virhe')
print -djpeg 'h6t4_virhe.jpg'
888888888888 88 8888888888
88 88 ,d 88
88 88 88 88 ____
88 ,adPPYba, 88,dPPYba, MM88MMM 88a8PPPP8b,
88 a8P_____88 88P' "8a 88 PP" `8b
88 8PP""""""" 88 88 88 d8
88 "8b, ,aa 88 88 88, 888 Y8a a8P
88 `"Ybbd8"' 88 88 "Y888 888 "Y88888P"
%tehdään koodi, joka tekee n. asteen approksimaatiot kaikilla tavoilla ja piirtää kuvan ja laskee virheet
>> x=linspace(-2,2);
>> a=-2;
>> b=2;
>> p=(x-a)/(b-a);
>> n=2;
F=zeros(1,n+1);
for k=0:n
F=F+h6t5(k/n)*bernstein(n,k,n+1);
end
B=F;
>> le=zeros(1,n+1);
for k=0:n
le=le+lege01(k,n+1)*kerroin(k);
end
>> k=n+1;
d=linspace(0,1,k);
>> d2=(d-a)/(b-a);
>> D=h6t5(d2);
la=polyfit(d2,D,n);
>> ph=a+(b-a)*x;
>>
x=linspace(0,1);
a=-2;
b=2;
x=phi(a,b,x);
n=2;
B=zeros(1,n+1); %bernstein B
for k=0:n
B=B+h6t5(phi(a,b,k/n))*bernstein(n,k,n+1);
end
le=zeros(1,n+1); %legendre le
for k=0:n
le=le+lege01(k,n+1)*kerroinh6t5(k);
end
k=n+1; %lagrange la
d=linspace(a,b,k);d2=linspace(0,1,k);
D=h6t5(d);
la=polyfit(d2,D,n);x=linspace(0,1);
yB=polyval(B,x); %arvoja
yle=polyval(le,x);
yla=polyval(la,x);
yf=h6t5(phi(a,b,x));
plot(x,yf) %kuvia
hold on
plot(x,yB,'r')
plot(x,yle,'k')
plot(x,yla,'m')
legend('f','B_n','le_n','la_n')
title('H6t5, f ja approksimaatiot')
vB=virheh6t5(B,10^3); %virheet
vle=virheh6t5(le,10^3);
vla=virheh6t5(la,10^3);
v=[vB,vle,vla]'
888888888888 88 ad8888ba,
88 88 ,d 8P' "Y8
88 88 88 d8
88 ,adPPYba, 88,dPPYba, MM88MMM 88,dd888bb,
88 a8P_____88 88P' "8a 88 88P' `8b
88 8PP""""""" 88 88 88 88 d8
88 "8b, ,aa 88 88 88, 888 88a a8P
88 `"Ybbd8"' 88 88 "Y888 888 "Y88888P"
%n. asteen tasavälinen interpolointipolynomi
n=10;
k=n+1;
x=linspace(-3,3,k);
y=h6t6(x);
la=polyfit(x,y,n);
X=linspace(-3,3);
Yf=h6t6(X);
Yla=polyval(la,X);
plot(X,Yf)
hold on
plot(X,Yla,'k')
eval(['legend(''f'',''tasavälinen ' num2str(n) '. asteen interpolointipolynomi'')']);
title('H6t6, f ja approksimaatio')
print -djpeg 'h6t6_tasavälinen.jpg'
%Rungen ilmiöstä voi tehdä videon seuraavasti.
for m=1:20 %piirretään m. asteisia polynomeja
n=m;
k=n+1;
x=linspace(-3,3,k);
y=h6t6(x);
la=polyfit(x,y,n);
X=linspace(-3,3);
Yf=h6t6(X);
Yla=polyval(la,X);
plot(X,Yf) %piirretään f
hold on
plot(X,Yla,'k') %piirretään polynomi
axis([-3 3 -3 3]) %laitetaan kuvaikkunan raamit aina näin
F(m)=getframe; %poimitaan ruutu
close %suljetaan kuvaikkuna, jotta saadaan puhdasta tilaa seuraavalla kierroksella
end
%n. asteen Tshebyshevin interpolointipolynomi
n=10;
for k=0:n
x(1+k)=3*cos(((2*k+1)*pi)/(2*n+2)); %Tshebyshevin pisteet
end
y=h6t6(x);
la=polyfit(x,y,n);
X=linspace(-3,3);
Yf=h6t6(X);
Yla=polyval(la,X);
plot(X,Yf)
hold on
plot(X,Yla,'k')
eval(['legend(''f'',''Tshebyshevin ' num2str(n) '. asteen interpolointipolynomi'')']);
title('H6t6, f ja approksimaatiot')
%kumpikin
n=10;
k=n+1; %tasavälinen
x=linspace(-3,3,k);
y=h6t6(x);
la=polyfit(x,y,n);
for k=0:n %Tshebyshev
x(1+k)=3*cos(((2*k+1)*pi)/(2*n+2)); %Tshebyshevin pisteet
end
y=h6t6(x);
ts=polyfit(x,y,n);
X=linspace(-3,3);
Yf=h6t6(X);
Yla=polyval(la,X);
Yts=polyval(ts,X);
plot(X,Yf)
hold on
plot(X,Yla,'k')
plot(X,Yts,'m')
eval(['legend(''f'',''tasavälinen ' num2str(n) '. asteen ip'',''Tshebyshevin ' num2str(n) '. asteen ip'')']);
title('H6t6, f ja approksimaatiot')