%esercizio 1 prova in itinere 26/01/2023
clear all
close all
clc

n=4;
x=linspace(0,2,8);
f=@(x) cos(x)-x.^2 -x +4;
y=f(x);

[C] = SISTC(x,n);
[Q,R]=QRC(C);
[a]=solvetriu(Q,R,y,n);

g=@(x) P4(x,a);
fplot(f,[-2,4])
hold on
fplot(g,[-2,4])

legend('Approssimazione','Funzione data')
t=linspace(3,4);
h=@(t) f(t) -P4(t,a) +1;
figure
fplot(h,[3,4])
yline(0,'--')

%circa tra 3.7 e 3.9

Nmax=100;
toll=sqrt(eps);
a=3.7;
b=3.9;

[t0,e,k]=METODOSECANTI(h,a,b,Nmax,toll);

hold on
plot(t0,0,'*r')

legend('f(x)-P4(x) +1','y=0','x_0');
xlabel('3 \leq x \leq 4')
ylabel('h(x)')

K=linspace(0,k,k+1);
figure
plot(K,e)
hold on
plot(K,e,'*r');

function [C] = SISTC(x,n)
%Costruisco C del mio sistema
m=length(x);
x=x';
C(:,[1 2])=[ones(m,1) x];

for i=1:m
    for j=3:n+1
        C(i,j)=C(i,j-1)*x(i);
    end
end
end

function [Q,R] = QRC(C)
%decomposizione QR
    [m,n]=size(C);
    Q=eye(m);

    for i=1:n
            v= C(i:m,i);
            u=zeros(m,1);
            if v(1)==0
                u(i:m)= v+ norm(v)*eye(m-i,1);
            else
                u(i:m)= v+ sign(v(1))*norm(v)*eye(m-i+1,1);
            end

            Qi=eye(m) -2*(u*u')/(u'*u);
            Q=Qi*Q;
            C=Qi*C;
    end    
    R=C;
end

function [a]=solvetriu(Q,R,y,n)
n=n+1;

b=Q*y';
b=b(1:n);
R=R(1:n,1:n);
a=eye(n,1);

a(n)=b(n)/R(n,n);

for i=n-1:-1:1
    a(i)=(b(i) - R(i,i+1:n)*a(i+1:n))/R(i,i);
end

end

function [P] = P4(xx,a)
%algoritmo valutazione polinomio 

n=length(a);
m=length(xx);

P=zeros(1,m);

for i=1:n
    P= P + a(i)*xx.^(i-1);
end
end

function [x0,e,k]=METODOSECANTI(h,a,b,Nmax,toll)
% metodo delle secanti
k=1;
x0=a;
x1=b;
r=(h(x1)-h(x0))/(x1-x0);
x1=x0;
x0=x0 - h(x0)/r;
e(1)=abs(h(x0));

while (k<Nmax) && (e(k)>toll)
    
    r=(h(x1)-h(x0))/(x1-x0);
    x1=x0;
    x0=x0 - h(x0)/r;
   
    k=k+1;
    e(k)=abs(h(x0));
end
k=k-1;
end