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 bernstein(n,k,N) antaan Bernsteinin polynomin B^n_k vähintään N pituisena vektorina. Esim.

bernstein(6,4,0)

ans =

    15   -30    15     0     0     0     0

bernstein(6,4,10)

ans =

     0     0     0    15   -30    15     0     0     0     0


%Nyt voidaan piirtää esimerkkitapaus. Vaikkapa

n=7;
k=3;
a=bernstein(n,k,0);
b=bernstein(n-1,k,0);
c=bernstein(n-1,k-1,0);
x=linspace(0,1);
A=polyval(a,x);
B=polyval(b,x);
C=polyval(c,x);
plot(x,A)
hold on
plot(x,B,'r')
plot(x,C,'k')
legend('B^5_4','B^4_4','B^4_3')
title('Harjoitus 5, tehtävä 2')

%Tämä voidaan tallentaa komennolla

print -djpeg 'h5t2b.jpg'


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'  


%Funktio f on tallennettu funktioksi h5t3.m



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')

%Pisteille saadaan kirjoitettua tehtävänannon mukaiset nimet komennolla

for k=0:n
eval(['text(K(' num2str(k+1) ',1),K(' num2str(k+1) ',2),''p' num2str(k) ''')']);
end


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"

f=[-3 4 5];	%tehtävänannon polynomi

L0=[0 0 1];	%kirjoitetaan Legendren polynomit L0-L2 niin, että ne esitetään
L1=[0 2 -1];	%samanpituisilla vektoreilla, jolloin niitä voi laskea yhteen.
L2=[6 -6 1];

x=linspace(0,1);	%valmistaudutaan integroimaan f*Li:tä välin 0..1 yli
fy=polyval(f,x);

y0=fy.*polyval(L0,x);		%f*L0
c0=(2*0+1)*trapz(x,y0);		%sen integraali kerrottuna kaavan vakiolla

y1=fy.*polyval(L1,x);
c1=(2*1+1)*trapz(x,y1);

y2=fy.*polyval(L2,x);
c2=(2*2+1)*trapz(x,y2);

p=c0*L0+c1*L1+c2*L2	%kirjoitetaan etsitty polynomi

p =

   -2.9847    3.9849    5.0024		%p on melkein f, joten kaikki meni hyvin

c=[c0 c1 c2]'				%etsityt vakiot ovat

c =

    5.9999
    0.5001
   -0.4974


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"   

%Tehdään 2. asteen approksimaatio Legendren polynomeilla kuten edellisessä tehtävässä
%nyt funktiona on kuitenkin h5t3.m

L0=[0 0 1];
L1=[0 2 -1];
L2=[6 -6 1];

x=linspace(0,1);
fy=h5t3(x);

y0=fy.*polyval(L0,x);		%f*L0
c0=(2*0+1)*trapz(x,y0);		%sen integraali kerrottuna kaavan vakiolla

y1=fy.*polyval(L1,x);
c1=(2*1+1)*trapz(x,y1);

y2=fy.*polyval(L2,x);
c2=(2*2+1)*trapz(x,y2);

p=c0*L0+c1*L1+c2*L2;		%etsitty polynomi
c=[c0 c1 c2]'			%etsityt kertoimet

py=polyval(p,x);

plot(x,fy)
hold on
plot(x,py,'r')

% verrataan tehtävän 3 approksimaatioon f2

n=2;
F=zeros(1,n+1);
for k=0:n
F=F+h5t3(k/n)*bernstein(n,k,n+1);
end
y=polyval(F,x);
plot(x,y,'k')

legend('h5t3 funktio f','approksimaatio Legendren polynomeilla','approksimaatio Bernsteinin polynomeilla')

c =

    1.2731
   -1.2160
   -1.3747