%=== nettoyage
clear all;
close all;
clc;

%=== parametres
n=16;%nombre de points a interpoler
nn=20;%nb de pts pour tracer chacune des splines

%=== points initiaux
x=sort(rand(n,1));%x croissants aleatoires entre 0 et 1
x(1,1)=0;
x(n,1)=1;
y=zeros(n,1);
for i=1:n
    xi=x(i,1);
    y(i,1)=cos(xi*8*pi)*xi;
end

%=== initialisation matrice vecteurs
a=zeros(4*(n-1),1);%[a1;b1;c1;d1;a2;b2;c2;d2;...]
b=zeros(4*(n-1),1);
P=zeros(4*(n-1),4*(n-1));

%=== remplissage matrice vecteur
for i=1:n-2
    xi=x(i+1,1);
    p=[1,xi,xi*xi,xi*xi*xi];
    %si(xi)=yi
    P(4*(i-1)+1,4*i+1:4*i+4)=p;
    b(4*(i-1)+1,1)=y(i+1,1);
    %si-1(xi)=yi
    P(4*(i-1)+2,4*(i-1)+1:4*(i-1)+4)=p;
    b(4*(i-1)+2,1)=y(i+1,1);
    
    p=[0,1,2*xi,3*xi*xi];
    %si'(xi)=si-1'(xi)
    P(4*(i-1)+3,4*i+1:4*i+4)=-p;
    P(4*(i-1)+3,4*(i-1)+1:4*(i-1)+4)=p;
    
    p=[0,0,2,6*xi];
    %si''(xi)=si-1''(xi)
    P(4*(i-1)+4,4*i+1:4*i+4)=-p;
    P(4*(i-1)+4,4*(i-1)+1:4*(i-1)+4)=p;
end

xi=x(1,1);
%sn-1(xn)=yn
p=[1,xi,xi*xi,xi*xi*xi];
P(4*(n-1)-3,1:4)=p;
b(4*(n-1)-3,1)=y(1,1);
%sn-1''(xn)=0
p=[0,0,2,6*xi];
P(4*(n-1)-2,1:4)=p;

xi=x(n,1);
%sn-1(xn)=yn
p=[1,xi,xi*xi,xi*xi*xi];
P(4*(n-1)-1,4*(n-1)-3:4*(n-1))=p;
b(4*(n-1)-1,1)=y(n,1);
%sn-1''(xn)=0
p=[0,0,2,6*xi];
P(4*(n-1)-0,4*(n-1)-3:4*(n-1))=p;

%=== calcul des coefficients des splines
a=P\b;

%=== trace
figure(1)
hold on;
plot(x,y,'ko','markerfacecolor',[1,1,1]);%points initiaux
for i=1:n-1
    xi=zeros(nn,1);
    yi=zeros(nn,1);
    yexact=zeros(nn,1);
    for j=1:nn
        xij=x(i,1)+(x(i+1,1)-x(i,1))/(nn-1)*(j-1);
        xi(j,1)=xij;
        for k=1:4
            yi(j,1)=yi(j,1)+a(4*(i-1)+k)*xij^(k-1);
        end
        yexact(j,1)=cos(xij*8*pi)*xij;
    end
    plot(xi,yexact,'-k','linewidth',2);
    plot(xi,yi,'-b','linewidth',2);
end

grid;
title('Interpolation 1D par splines cubiques f(x)=x*cos(8*pi*x)');
legend('points','exact','spline','location','northwest');
xlim([0,1]);
hold off;