clear;
%--------------------------------------
% Esercizio 1 del 6/3/06
% Approssimazione della derivata prima
% di una funzione
%--------------------------------------
a=0;
b=pi;
n=input('Numero dei punti della suddivisione = ');
tipoSchema=input('Quale schema? A=DiffFinAvanti, I=DiffFinIndietro e C=DiffFinCentrali  ','s');
index=input('In quale punto INTERNO desideri valutare la derivata prima?  ');
h=(b-a)/n;
x=a:h:b;
y=funD(x);

xx=a:0.01:b;
[yy,yy1]=funD(xx);
figure(1);
plot(xx,yy,'-',x,y,'or')
legend('Funzione', 'Punti equispaziati');
figure(2)
plot(xx,yy1,'-',x,y,'or')
legend('Derivata funzione', 'Punti equispaziati');

fpA=funD(x(index+1));
[fAi,fDer]=funD(x(index)); %qui calcolo anche il valore della derivata VERA
fi=funD(x(index-1));

xx=x(index):0.001:x(index+1);
[A,B,f2A,C]=funD(xx);
xx=x(index-1):0.001:x(index);
[A,B,f2I,C]=funD(xx);
xx=x(index-1):0.001:x(index+1);
[A,B,C,f3]=funD(xx);

switch tipoSchema
    case {'A'}
        disp('Derivata approssimata ')
        fdA=(fpA-fAi)/h
        disp('Derivata vera ');
        fDer
        disp('Stima dell'' errore di approssimazione')
        errA=h^2/2*max(abs(f2A))
    case {'I'}
        disp('Derivata approssimata ')
        fdI=(fAi-fi)/h
        disp('Derivata vera ');
        fDer
        disp('Stima dell'' errore di approssimazione')
        errI=h^2/2*max(abs(f2I))
    case {'C'}
        disp('Derivata approssimata ')
        fdC=(fpA-fi)/(2*h)
        disp('Derivata vera ');
        fDer
        disp('Stima dell'' errore di approssimazione')
        errC=h^2/6*max(abs(f3))
    otherwise 
        error('Tipo schema non previsto o errata scrittura')
end

function [y,yD,yDD,yDDD]=funD(x)
sx=sin(x.^2)
cx=cos(x.^2);
y=cos(sx);
yD=-2*x.*sin(sx).*cx;
yDD=-4*cos(sx).*cx.^2.*x.^2+4*sin(sx).*sx.*x.^2-2*sin(sx).*cx;
yDDD=4*x.*(2*sin(sx).*cx.^3.*x.^2+6*cos(sx).*cx.*x.^2.*sx-3*cos(sx).*cx.^2+2*sin(sx).*cx.*x.^2+3*sin(sx).*sx);
return

%-------------------------------------------------------------------------------------

clear;
%---------------------------------
% Esercizio 2 del 06/03/06 su
% integrazione numerica
%---------------------------------
format long
tipoFormula=input('Quale formula: (T) trapezi, (S) Simpson o (G) gaussiana ? ','s');
n=input('Numero di punti della suddivisione ');
a=0;
b=2*pi;
h=(b-a)/n;
x=linspace(a,b,n+1); %punti equispaziati.
xx=a:0.01:b;
yy=funQ(xx);
realValue=quadl(@funQ,a,b,1.e-6);
switch tipoFormula
    case {'T'}
        fTc=funQ(x);
        plot(xx,yy,'-.g',x,fTc,'o');
        fTc(2:end-1)=2*fTc(2:end-1);
        ValTc=0.5*h*sum(fTc)
        title('Quadratura composita con i trapezi e relativi nodi');
        disp('Errore assoluto')
        erroreT=abs(realValue-ValTc)
    case {'S'}
        fSc=funQ(x);
        plot(xx,yy,'-.g',x,fSc,'ob'), hold on
        fSc(2:end-1)=2*fSc(2:end-1);
        ValSc=h*sum(fSc)/6;
        x=linspace(a+h/2,b-h/2,n); 
        fSc=funQ(x);
        plot(x,fSc,'or')
        ValSc=ValSc+2*h/3*sum(fSc)
        title('Quadratura composita di Simpson e relativi nodi');
        disp('Errore assoluto')
        erroreT=abs(realValue-ValSc)
        hold off
    case {'G'}
        y1=x(1:end-1)+h/2*(1-1/sqrt(3));
        plot(xx,yy,'-.g',y1,funQ(y1),'ok'); hold on
        y=[y1 y1+h/sqrt(3)];
        fGc=funQ(y);
        plot(xx,yy,'-.g',y,fGc,'or')
        ValGc=0.5*h*sum(fGc)
        title('Quadratura composita con formula di Gauss e relativi nodi');
        disp('Errore assoluto')
        erroreT=abs(realValue-ValGc)
        hold off
    otherwise
        error('Formula sconosciuta')
end

function y=funQ(x)
y=x.*exp(-x).*cos(2*x);
end
