%es 1 23.07.21

clear all
close all
clc

x=[-3 -2 -1 0 1 2];
y=[5 2.5 0.5 0 -0.4 -0.7];
xx=linspace(x(1),x(end));

%PUNTO A
% Sarà un polinomio di grado 5, utilizzo Lagrange

P=@(xx) Lagrange(x,y,xx);

%PUNTO B
figure
subplot(2,1,1)
fplot(P)
title('Approssimazione Polinomiale','Su R')
hold on
subplot(2,1,2)
hold on
fplot(P,[x(1),x(end)])
plot (x,y,'*')
legend('Metodo di Lagrange','(x,f(x)')
xlabel('-3 leq x \leq 2')
ylabel ('y(x)')
title(' ','Zoom su intervallo di interesse')

%PUNTO C

f=@(xx) P(xx)-exp(xx)+1;

%Newton-Raphoson mi assicurerebbe una velocità di convergenza maggiore
% dovrei però vedere se si annulla la derivata  nell'intervallo di ricerca,
% uso il metodo Newton alle differenze e confronto la convergenza col
% metodo delle secanti (velocità elevata)

Nmax=100;
toll=sqrt(eps);

[xd,ed] = DiffFin(f,xx,Nmax,toll);
[xs,es] = Secanti(f,xx,Nmax,toll);

figure
fplot(f,[xx(1),xx(end)])
hold on
plot(xd,f(xd),'*b')
plot(xs,f(xs),'*r')
title('Ricerca 0')

figure
subplot(2,1,1)
plot(ed,'g')
hold on
plot(es,'b')
title('Convergenza Errore','Forma canonica','FontAngle','italic')

subplot(2,1,2)
plot(log10(ed),'g')
hold on
plot(log10(es),'b')
title('','Forma logaritmica','FontAngle','italic')

function [yy] = Lagrange(x,y,xx)
n=length(x)-1;
yy=xx-xx;

for j=0:n
    num=1;
    den=1;
    for i=[0:j-1,j+1:n]
        num=num.*(xx-x(i+1));
        den=den.*(x(j+1)-x(i+1));
    end
    yy=yy + y(j+1).*num./den;
end

end
function [xd,ed] = DiffFin(f,xx,Nmax,toll)

xd=xx(1);
if abs(xd)<toll
        h=toll;
    else
        h=toll*xd;
    end
r=(f(xd+h)-f(xd))/h;
xd= xd - f(xd)/r;
k=1;
ed(k)=abs(f(xd));

while k<Nmax && ed(k)>toll
    if abs(xd)<toll
        h=toll;
    else
        h=toll*xd;
    end
    r=(f(xd+h)-f(xd))/h;
    xd= xd - f(xd)/r;
    k=k+1;
    ed(k)=abs(f(xd));
end

end

function [xs,es] = Secanti(f,xx,Nmax,toll)
xs=xx(1);
x=xx(end);
r=(f(xs)-f(x))/(xs-x);
x=xs;
xs= xs -f(xs)/r;
k=1;
es(k)=abs(f(xs));

while k<Nmax && es(k)>toll
    r=(f(xs)-f(x))/(xs-x);
    x=xs;
    xs= xs -f(xs)/r;
    k=k+1;
    es(k)=abs(f(xs));
end

end