%% Harj. 2, teht. 6 ja meshgrid-selvitysta
%
%% meshgrid, ulkotulo 
% Liittyy harj.2 teht. 6,7
% Aja t"am"a kappale, niin ymm"arr"at meshgridin.
% Hy"odynnet"a"an Matlabin symbolic toolboxia.
clear
close all
syms x1 x2 x3
syms y1 y2 y3 y4 y5
x=[x1 x2 x3]
y=[y1 y2 y3 y4 y5]

[X,Y]=meshgrid(x,y)

Z=X.*Y

%%

% Sama syntyisi for-for-silmukalla:
% (Mutta olispa tyhm"a"a, jos pit"ais rakentaa t"allainen aina,
%  kun haluaisi 3d-grafiikkaa tuottaa.)
clear W
for i=1:5
    for j=1:3
        W(i,j)=x(j)*y(i);
    end
end    
W
%%
% Korkeusmatriisi funktiolle f(x,y)=cos(x^2)sin(y^2)

Z1=cos(X.^2).*sin(Y.^2)

%% Teht. 7
%
xd=linspace(0,3,7);
yd=cos(1+xd.^2);

c=polyfit(xd,yd,6)   % 7 pistett"a => interpolaatiopolynomin asteluku=6
x=linspace(0,3);     % Pisteet, joissa interpolaatiopoly halutaan laskea
p=polyval(c,x);      % No, lasketaan.
y=cos(1+x.^2);
plot(x,p,'b',x,y,'r',xd,yd,'ok') % Piirret"a"an polynomi,funktio ja datapisteet
grid on
legend('Interpolaatiopolynomi','funktio','data','Location','SouthWest')
shg
% Lasketaan samantien todellinen max-virhe (100:n pisteen resoluutiolla):
TodellinenVirhe = max(abs(y-p))
%%

% Selvitet"a"an ensin alla oleva periaate.

%clear 
close all
syms xd0 xd1 xd2
syms x1 x2 x3 x4 x5
xd=[xd0 xd1 xd2]
x=[x1 x2 x3 x4 x5]

[X,Xd]=meshgrid(x,xd)
X-Xd
prod(X-Xd)
abs(prod(X-Xd))

%% Virhetermi: ensin tulotermin max
% Edell"a malli.
xd=linspace(0,3,7);
yd=cos(1+xd.^2);
x=linspace(0,3);
[X,Xd]=meshgrid(x,xd);
tulot=prod(X-Xd);
maxtulot100=max(abs(tulot)) % Resoluutiolla 100 pistett"a.
%{ 
maxtulot100 =
    0.7470
Vrt: Maplella analyyttisesti: 0.7487
%}
plot(x,abs(tulot));shg
%
%%
% Kokeillaan huvin vuoksi resoluutiolla 200:
x=linspace(0,3,200);
[X,Xd]=meshgrid(x,xd);
tulot=prod(X-Xd);
maxtulot200=max(abs(tulot)) % Resoluutiolla 200 pistett"a.
%{
maxtulot200 =
    0.7483
Saatiin 3 oikeaa numeroa. 
T"am"an teht"av"an kannalta ei merkityst"a, mutta yleisesti kiintoisaa.
Toisinaan "raakaan voimaan" perustuva laskenta riitt"a"a hienon analyysin
sijasta, t"ass"a nyt ilman muuta.
%}
maxT=maxtulot200
%% 7. derivaatta
syms x
d7=diff(cos(1+x^2),x,7)
%{
d7 =
1680*x*cos(x^2 + 1) - 1344*x^5*cos(x^2 + 1) - 3360*x^3*sin(x^2 + 1) + 128*x^7*sin(x^2 + 1)
%}
d7=vectorize(d7)  % Laskutoimitukset vektoriop:ksi, eli lis"at"a"an pisteet.
%{
d7 =
1680.*x.*cos(x.^2 + 1) - 1344.*x.^5.*cos(x.^2 + 1) - 3360.*x.^3.*sin(x.^2 + 1) + 128.*x.^7.*sin(x.^2 + 1)
%}
x=linspace(0,3);
d7num=eval(d7);  % Huom! d7 ei p"aivity itsest"a"an, kuten Maplessa
plot(x,d7num)
maxd7=max(abs(d7num))
%{
maxd7 =
  2.7618e+005
Maplesta poimittu: 2.76*10^5
Hyvin "matchaa"!
%}
%%   Virhearvio maksimivirheelle:
maxvirhearvio=maxd7/prod(1:7)*maxT
%{
maxvirhearvio =
   41.0092
Kuten Maplessa
%}
%% Todellinen virhe
%Laskettiin yll"a:
TodellinenVirhe
%{
TodellinenVirhe =
    0.4535

N"ainkin huonosti voi menn"a, virhearvio on 100-kertainen todelliseen
verrattuna. No, ongelma on tuo tolkuton 7. derivaatan maksimi.
%}
