// 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 résoud un système d'équation
// pour déterminer les polynômes d'interpolation.

// 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

// Taille par défaut des fenêtres.
windef = gdf();
windef.figure_size = [580, 420];
windef.axes_size = [550,400];

function y=fct(x)
//===============
// La fonction dont on a pris des points d'interpolations.
y = 0.5 + erf(x /sqrt(2))/2;
endfunction

ax = linspace(0, 2, 6);  // 6 points d'interpolation
ay = [0.5, 0.6554217416, 0.7881446014, 0.8849303298, 0.9452007083, 0.9772498681];
//ay = fct(ax) + 2e-3 * (rand(ax) - 0.5);  // on ajoute du bruit sur les données !

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 y=poly_f(coef, x)
//========================
// Polynome de degré  length(coef)-1,  de coefficients coef
pp = ones(x) * coef($);
for ii=length(coef)-1:-1:1 
   pp = pp .* x + coef(ii);
   end;
y = pp;
endfunction

function z=interpol_poly(axn, ayn, x)
//===================================
// Interpolation du polynôme de degré length(axn)-1
// passant pas les points  axn,  ayn
// calcul des coefficients du polynôme passant par ax2,  ay2

// Il n'est pas efficace de résoudre le système à chaque fois
// mais le but n'est pas d'être efficace ici, mais de se simplifier la vie.

// Système d'équations à résoudre.
A = ones(axn)';
for ii=1:length(axn)-1;  A = [A, (axn.^ii)']; end;

// Coefficient du polynome passant par les points défini par axn, ayn
coef = A\ayn';

z = poly_f(coef, x);
endfunction

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

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

// Graphique de l'erreur de l'interpolation polynomiale de degré  length(ax)-1
scf(5);  // 5 est le degré du polynôme dans cet exemple
clf();
plot(xxx, interpol_poly(ax, ay, 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 de degré 2
//==========================
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) = interpol_poly(ax(j-1:j+1), ay(j-1:j+1), x(i));   
end; // for
endfunction

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

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

// Interpolation de degré 3
//==========================
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) = interpol_poly(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 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




