// Solution de la série 3 sur l'interpolation
// L'objectif est d'interpoler une fonction, connue que en quelques points

// Dans cet exemple, on utilise les différences divisées pour évaluer des polynôme.

// Référence : http://www.unige.ch/~hairer/poly/poly.pdf  
// Chapitre II, Interpolation et Approximation
// Plus spécifiquement : chapitres II.1, II.2, (II.3), II.10

funcprot(0); // pour ne pas indiquer une redéfinition des fonctions

windef = gdf();
windef.figure_size = [580, 420];
windef.axes_size = [550,400];

function y=fct(x)
//===============
y = 0.5 + erf(x /sqrt(2))/2;
endfunction

ax = linspace(0, 2, 6);
ay = [0.5, 0.6554217416, 0.7881446014, 0.8849303298, 0.9452007083, 0.9772498681];

xxx=0:0.01:max(ax);  // pour les plot

function z=intlin(x1, x2, y1, y2, x)
//==================================
// Interpolation linéaire 
z = y1 + (y2 - y1) / (x2 - x1) * (x - x1);
endfunction

scf(0);
clf(0);

//plot(ax, ay, 'LineStyle', 'none', 'Marker', 'o');
plot(ax, ay, 'Marker', 'o');

// Définition de l'échelle des axes
mygca=gca();
mygca.zoom_box = [min(ax)-0.1, min(ay)-0.05, max(ax)+0.1, max(ay)+0.05];

x_milieu = (ax(3) + ax(4)) / 2; 
z1 = intlin(ax(3), ax(4), ay(3), ay(4), x_milieu);

// Le point d'interpolation au milieu du 3ème segment
plot(x_milieu, z1, 'Marker', 'o');
mypts = mygca.children(1).children; // structure de paramètres des points
mypts.mark_size=9; // taille des points

// Divers points d'interpolation linéaire du 3ème segment
ax2=ax(3):0.05:ax(4);
plot(ax2, intlin(ax(3), ax(4), ay(3), ay(4), ax2), 'Marker', 'o');
mypts = mygca.children(1).children; // structure de paramètres des points
mypts.mark_size=5; // taille des points
mypts.mark_style=3;
title("Données");

//####################################################################

function fx=lineaire_interpol(ax, ay, x)
//======================================
// Retourne l'interpolation linéaire en  x  de la fonction
// passant par les points  ax, ay
fx = zeros(x);

for i=1:length(x)
  // recherche dans quel intervalle  [ax(j), ax(j+1)]  se trouve  x(i)
  j = length(ax) - 1;
  while (ax(j) > x(i)) & ( j > 1); j = j-1; end;

  // Ici  ax(j) <= x(i) < ax(j+1),  sauf si  x(i) > max(ax) ou x(i) < min(ax)
  fx(i) = intlin(ax(j), ax(j+1),ay(j), ay(j+1), x(i));   
end; // for
endfunction

// Graphique de l'erreur de l'interpolation linéaire
scf(1);
clf(1);
plot(xxx, lineaire_interpol(ax, ay, xxx) - fct(xxx));
title("Interpolation linéaire");

//###########################################################################

function az=differences_divisees(ax, ay)
//======================================
// Pour deux vecteurs  ax  et  ay,  
// Retourne le vecteur des différences divisées de f(ax) = ay
// C'est utile pour calculer un polynôme d'interpolation
az = ay;
for i = 2:length(ax)
  for j=length(ax):-1:i
    az(j) = (az(j) - az(j-1)) / (ax(j)- ax(j+1-i));
  end
end
endfunction

// Interpolation de degré 2
//==========================
function z=intquadratique(ax2, ay2, x)
//====================================
// Interpolation quadratique

// calcul des différences divisées de  ax2,  ay2
dd2 = differences_divisees(ax2, ay2);

//z = dd2(1) + (x - ax2(1)).*(dd2(2) + (x - ax2(2)).*dd2(3));
// Mieux
z = (dd2(3).*(x - ax2(2)) + dd2(2)).*(x - ax2(1)) + dd2(1);
endfunction

function fx=quadratique_interpol(ax, ay, x)
//=========================================
// Retourne l'interpolation quadratiquecubique en  x  de la fonction
// passant par les points  ax, ay
fx = zeros(x);

for i=1:length(x)
  // recherche dans quel intervalle  [ax(j), ax(j+1)]  se trouve  x(i)
  j = length(ax) - 1;
  while (ax(j) > x(i)) & ( j > 2); j = j-1; end;

  // Ici  ax(j) <= x(i) < ax(j+1),  sauf si  x(i) > max(ax) ou x(i) < min(ax)
  fx(i) = intquadratique(ax(j-1:j+1), ay(j-1:j+1), x(i));   
end; // for
endfunction

// Graphique de l'erreur de l'interpolation cubique
scf(2);
clf();
int2 = quadratique_interpol(ax, ay, xxx);
plot(xxx, int2 - fct(xxx));
title("Interpolation quadratique");

//###########################################################################

// Interpolation de degré 3
//==========================
function z=intcubique(ax3, ay3, x)
//================================
// Interpolation cubique

// calcul des différences divisées de  ax3,  ay3
dd3 = differences_divisees(ax3, ay3);

//z = dd3(1) + (x - ax3(1)).*(dd3(2) + (x - ax3(2)).*(dd3(3) + (x - ax3(3))*dd3(4)));
// Mieux
z = ((dd3(4).*(x - ax3(3)) + dd3(3)).*(x - ax3(2)) + dd3(2)).*(x - ax3(1)) + dd3(1);
endfunction

function fx=cubique_interpol(ax, ay, x)
//=====================================
// Retourne l'interpolation cubique en  x  de la fonction
// passant par les points  ax, ay
fx = zeros(x);

for i=1:length(x)
  // recherche dans quel intervalle  [ax(j), ax(j+1)]  se trouve  x(i)
  j = length(ax) - 2;
  while (ax(j) > x(i)) & ( j > 2); j = j-1; end;

  // Ici  ax(j) <= x(i) < ax(j+1),  sauf si  x(i) > max(ax) ou x(i) < min(ax)
  fx(i) = intcubique(ax(j-1:j+2), ay(j-1:j+2), x(i));   
end; // for
endfunction

// Graphique de l'erreur de l'interpolation cubique
scf(3);
clf();
int3 = cubique_interpol(ax, ay, xxx);
plot(xxx, int3 - fct(xxx));
title("Interpolation cubique");

//###########################################################################

// Interpolation de degré 5

function z=intcinq(ax5, ay5, x)
//================================
// Interpolation de degré 5

// calcul des différences divisées de  ax5,  ay5
dd5 = differences_divisees(ax5, ay5);

//z = dd5(1) + (x - ax5(1)).*(dd5(2) + (x - ax5(2)).*(dd5(3) + ..
//    (x - ax5(3))*(dd5(4) + (x - ax5(4)).*(dd5(5) + (x - ax5(5)).*dd5(6)))));
// MIEUX
z = ((((dd5(6).*(x - ax5(5)) + dd5(5)).*(x - ax5(4)) + dd5(4)).* ..
    (x - ax5(3)) + dd5(3)).*(x - ax5(2)) + dd5(2)).*(x - ax5(1)) + dd5(1);
    
endfunction

function fx=cinq_interpol(ax, ay, x)
//==================================
// Retourne l'interpolation de degré 5 en  x  de la fonction
// passant par les points  ax, ay
fx = zeros(x);

for i=1:length(x)
  // recherche dans quel intervalle  [ax(j), ax(j+1)]  se trouve  x(i)
  j = length(ax) - 3;
  while (ax(j) > x(i)) & ( j > 3); j = j-1; end;

  // Ici  ax(j) <= x(i) < ax(j+1),  sauf si  x(i) > max(ax) ou x(i) < min(ax)
  fx(i) = intcinq(ax(j-2:j+3), ay(j-2:j+3), x(i));   
end; // for
endfunction

// Graphique de l'erreur de l'interpolation de degré 5
scf(5);
clf();
int5 = cinq_interpol(ax, ay, xxx);
plot(xxx, int5 - fct(xxx));
title("Interpolation de degré 5");

//###########################################################################

// calcul des différences divisées de  ax,  ay
difd = differences_divisees(ax, ay);

function fx=polynome_interpol(ax, ad_d, x)
//========================================
// Pour  ax  et les différences divisées  ad_d,  
// évalue le polynome d'interpolation en  x,  x pouvant être un vecteur
// Le retour  fx  est un vecteur de même dimension que  x
//fx = ad_d(1) * ones(x);
//terme = ones(x);
//for i=2:length(ax);
//  terme = terme .* (x - ax(i-1)); 
//  fx = fx + terme * ad_d(i);
//end

// MIEUX
fx = ad_d(length(ax)) * ones(x);
for i=length(ax)-1:-1:1;
  fx = fx .* (x - ax(i)) + ad_d(i);
end

endfunction

// Montre que le polynome d'interpolatoin pass bien par les points ax, ay
if %F then
  for i=1:length(ax)
    f1 = polynome_interpol(ax, difd, ax(i));
    disp(f1 - ay(i));
  end; // for
end; // if

// Graphique de l'erreur de l'interpolation polynomiale de degré  length(ax)-1
scf(6);
clf();
plot(xxx, polynome_interpol(ax, difd, xxx) - fct(xxx));
title("Interpolation polynomiale de degré " + string(length(ax)-1));
//xstring(0, 0, "l''interpolation polynomiale est bonne jusqu''au degré 9");
xlabel("L''interpolation polynomiale est bonne jusqu''au degré 9");

//###########################################################################

// Interpolation spline
//=====================
// Cas du spline naturel, qui est assez bon
derivees = splin(ax, ay, "natural");

// Graphique de l'erreur de l'interpolation spline
scf(7);
clf();
spline1 = interp(xxx, ax, ay, derivees);
plot(xxx, spline1 - fct(xxx));
title("Interpolation spline");

// Référence : http://www.unige.ch/~hairer/poly/poly.pdf  
// Chapitre II, Interpolation et Approximation
// Chapitre II.10, Interpolation par fonction spline




