%esercizio 2 19.09.2022
close all
clear all
clc

x=linspace(0,3,10)';
xx=linspace(x(1),x(end));
f=@(xx) -6*sin(2*xx) + xx.^2 +7*xx+1;
y=f(x);

g=@(xx) spline(x,y,xx);
yy=spline(x,y,xx);
fplot(g,[x(1),x(end)])
fplot(f,[x(1),x(end)])
hold on
plot (x,y,'*r')
plot(xx,yy,'*b')

% versione vecchia --> discontinuità
% 
% function [yy]=spline(x,y,xx)
% h=diff(x);
% v=diff(y)./h;
% 
% D=2*(x(3:end)-x(1:end-2)); 
% C=h(2:end-1);              
% 
% A=diag(D)+diag(C,-1)+diag(C,1);
% b=6*(v(2:end)-v(1:end-1));
% z=A\b; %risolvere il sistema sarebbe più bello!!
% z=[0; z; 0];
% 
% yy=zeros(size(xx));
% for ii=1:length(x)-1
%     I=(x(ii)<=xx)&(x(ii+1)>=xx);
%     
%     C=y(ii+1)./h(ii)-h(ii)/6*z(ii+1);
%     
%     D=y(ii)/h(ii)-h(ii)/6*z(ii);
%     
%     yy=yy+(z(ii+1)/6/h(ii).*(xx-x(ii)).^3+z(ii)/6/h(ii).*(x(ii+1)-xx).^3+...
%         +C*(xx-x(ii))+D*(x(ii+1)-xx)).*I;
% end
% end

% versione nuova, la vecchia creava delle discontinuità in prossimità dei
% miei nodi;
% La supposizione è che fosse dovuto alla doppia valutazione delle
% condizioni logiche del vettore I nei miei nodi interni.
% In prima linea ho cambiato quello, ma togliere da uno dei due estremi
% degli intervallini la possibile uguaglianza implica che io stia
% trascurando un nodo (minore o maggiore) esterno (a seconda di quale
% uguaglianza trascuro).
% Per ovviare a questo aggiungo un controllo di condizione per il mio
% indice. In questo caso alla prima iterazione posiziono un 1 nel mio
% vettore logico in prima posizione.
% Ci saranno sicuramente metodi più efficienti e forse non era proprio
% dovuto a quello --> chiedi Prof.ssa Martini 
% 
function [yy]=spline(x,y,xx)
h=diff(x);
v=diff(y)./h;

D=2*(x(3:end)-x(1:end-2));
C=h(2:end-1);              

A=diag(D)+diag(C,-1)+diag(C,1);
b=6*diff(v);
M=A\b;
M=[0; M; 0];

yy=zeros(size(xx));
for ii=1:length(x)-1

    if ii==1
        I=(x(ii)<=xx)&(x(ii+1)>=xx);
    else
        I=(x(ii)<xx)&(x(ii+1)>=xx);
    end

    A=M(ii+1)/6/h(ii).*(xx-x(ii)).^3;
    B=M(ii)/6/h(ii).*(x(ii+1)-xx).^3;
    C=y(ii+1)./h(ii)-h(ii)/6*M(ii+1);
    D=y(ii)/h(ii)-h(ii)/6*M(ii);
    
    yy=yy+(A+B+C*(xx-x(ii))+D*(x(ii+1)-xx)).*I;
end

end 

% % In questo caso invece lo metto in ultima posizione nell'ultima
% iterazione; 

% 
% function [yy]=spline(xx,x,y)
% h=diff(x);
% v=diff(y)./h;
% 
% D=2*(x(3:end)-x(1:end-2)); %costruisco D
% C=h(2:end-1);              %costruisco C
% 
% %Az=b
% A=diag(D)+diag(C,-1)+diag(C,1);
% b=6*diff(v);
% z=A\b;
% z=[0; z; 0];
% 
% %Costruisco yy
% yy=zeros(size(xx));
% for ii=1:length(x)-1
% 
%     if ii==length(x)-1
%         I=(x(ii)<=xx)&(x(ii+1)>=xx);
%     else
%         I=(x(ii)<=xx)&(x(ii+1)>xx);
%     end
% 
%     C=y(ii+1)./h(ii)-h(ii)/6*z(ii+1);
%     
%     D=y(ii)/h(ii)-h(ii)/6*z(ii);
%     
%     yy=yy+(z(ii+1)/6/h(ii).*(xx-x(ii)).^3+z(ii)/6/h(ii).*(x(ii+1)-xx).^3+...
%         +C*(xx-x(ii))+D*(x(ii+1)-xx)).*I;
% end
% 
% end