// Le but de ce qui suit est d'obtenir le "meilleur" polynome de degré donné
// approximant une fonction donnée sur un intervalle donné.
// Meilleure est dans le sens que :
// max(abs(f(x) - p(x)) = minimum,  
// a < x < b, f  étant la fonction et  p  le polynome
// L'algorithme de Remez permet de s'approcher rapidement de la meilleure
// approximation polynomiale.
// c.f. : https://fr.wikipedia.org/wiki/Algorithme_de_Remez
// c.f. : https://en.wikipedia.org/wiki/Remez_algorithm

// EN CAS DE PROBLEMES, taper  clear  dans la console.

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

// Définit les dimensions des fenêtres graphique
windef = gdf();
windef.figure_size = [580, 420];
windef.axes_size = [550,400];

function y=fct(x)
//===============
// La fonction à approximer
// Les approximations sont bonne pour les degrés allant de 5 à 10,
// Pour  nnn > 10, les limites de précision numérique sont parfois atteinte.
//y = 1 - erf(x * sqrt(2))/2; 
//y = 0.5 + erf(x * sqrt(2))/2;
//y = exp(-5*x.^2);
y = sin(x*%pi/2); // pour nnn > 10, les limites de précision sont atteintes
//y = sin(sqrt(x)*%pi/2) ./ sqrt(x);
//y = exp(x); // excellante approximation de degrés 4, 6, 8.  nnn=10 atteint les limites de précision
//y = cos(x*%pi/2); // pour nnn > 10, les limites de précision sont atteintes
//y = sqrt(1 + 3*x); // 
//y = 1 ./ sqrt(1 + 3*x); 
//y = 1 ./ (1 + 25*x.^2); // Dans aucun cas l'approximation est excellante !

//y = tan(x*%pi/2) .* (1 - x);
//y = x ./ tan(x*%pi/2);
  // La précision de cette approximation est la même que celle de la fonction qui précède
endfunction

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

//########### DONNEES #################

aaa = 1e-9;  // borne inférieure.  > 0 pour éviter des divisions par 0, selon la fonction
bbb = 1-1e-9;  // borne supérieure
nnn = 6; // nnn+1 = nombre de points   nnn = degré du polynôme approximant
xxx = linspace(aaa, bbb, 20*nnn+1);  // pour les plot

scf(0);
clf();
yyy = fct(xxx);
plot(xxx, yyy);

// Définition de l'échelle des axes
mygca=gca();
mygca.zoom_box = [aaa-0.1, min(yyy)-0.05, bbb+0.1, max(yyy)+0.05];
title("Courbe à approximer.");

//############ 1111111111111111111111111111 #################

// Détermination du polynôme d'interpolation passant par les points régulièrement espacés
kkk = 0:nnn;
ax = linspace(aaa, bbb, nnn+1); 
ay = fct(ax);

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

// Coefficient du polynome passant par les nnn+1 points défini par ax, ay
coef = A\ay';

scf(1);
clf();
plot(xxx, poly_f(coef, xxx) - fct(xxx));
title("Points régulièrement espacés. Erreur de l''interpolation polynomiale de degré " +..
      string(nnn));

//############ 2222222222222222222222222222 #################

// Détermination du polynôme d'interpolation passant par les points de Chebyshev (Tchebyshev)
// du polynôme de Tchebyshev de degré  nnn+1
kkk = 0:nnn;
ax = (aaa+bbb)/2 + (bbb-aaa)/2 * cos((2*kkk+1)*%pi / (2*nnn+2));
// Les points du vecteurs  ax  se nomment les points de Chebyshev
// Ils sont bons pour être proche de l'approximation désirée.
// c.f. "Polynômes de Chebyshev" du cours : http://www.unige.ch/~hairer/poly/poly.pdf
ay = fct(ax);

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

// Coefficient du polynôme passant par les nnn+1 points défini par ax, ay
coef = A\ay';

scf(2);
clf();
plot(xxx, poly_f(coef, xxx) - fct(xxx));
title("Aux Zéros du polynôme de Chebyshev. Erreur de l''interpolation polynomiale de degré " +..
      string(nnn));

//############ 3333333333333333333333333333 #################

// Meilleures points pour une approcher une approximation optimale
// Les points  xx  ci-dessous correspondent aux extremums de la différence
// entre l'approximation polynomiale et la fonction.
// Extreumus du polynôme de Tchebyshev de degré  nnn+1
kkk = 0:nnn+1;
xx = (aaa+bbb)/2 + (bbb-aaa)/2 * cos((kkk*%pi / (nnn+1)));
xx = xx'; // pour avoir un vecteur colonne

ay = fct(xx);
A = ones(ay);
for ii=1:nnn;  A = [A, xx.^ii]; end;
for ii=1:nnn+2;  A(ii, nnn+2) = (-1)^ii; end;

// Coefficient du polynôme passant proche de l'optimum
coef = A\ay;
disp("Aux extremums de Chebyshev, erreur max env. = " + string(coef($)));
coef = coef(1:$-1);

scf(3);
clf();
plot(xxx, poly_f(coef, xxx) - fct(xxx));
title("Aux extremums du polynôme de Chebyshev. Erreur de l''interpolation polynomiale de degré " +..
       string(nnn));

//############ 4444444444444444444444444444 #################

// Implémentation de l'algorithme de Remez
// =======================================

function [xm, ym]=xMinimum(coef, xExtr, yErr1)
//============================================
// Recherche du minimum de la fonction différence : dif(x) = poly_f(x) - fct(x)
// à partir de xExtr.  yErr1 = diff(xExtr)
// (xm, ym) = coordonnées du minimum 
yErr2 = 1;
while (yErr1 < yErr2) & (xExtr < bbb)
  yErr2 = yErr1;
  xExtr = xExtr + 0.001;
  yErr1 = poly_f(coef, xExtr) - fct(xExtr);
  end;
xm = xExtr;
ym = yErr1;  
endfunction

function [xm, ym]=xMaximum(coef, xExtr, yErr1)
//============================================
// Recherche du maximum de la fonction différence : diff(x) = poly_f(x) - fct(x)
// à partir de xExtr.  yErr1 = diff(xExtr)
// (xm, ym) = coordonnées du maximum
yErr2 = -1;
while (yErr1 >= yErr2) & (xExtr < bbb)
  yErr2 = yErr1;
  xExtr = xExtr + 0.001;
  yErr1 = poly_f(coef, xExtr) - fct(xExtr);
  end;
xm = xExtr;
ym = yErr1;  
endfunction

// Une ou deux itérations de l'algorithme de Remez suffisent
// nnn = degré du polynôme d'interpolation.
// nnn+2 = nombre d'extremas, y compris les deux extrémités de l'intervalle
for iter=1:1
  // Recherche des maximums, en 3 itérations
  xx(1) = aaa;
  xx(nnn+2) = bbb;
  xExtr = aaa;
  yErr1 = poly_f(coef, xExtr) - fct(xExtr);

  jj=1;
  while (jj <= nnn)
    if (yErr1 >= 0) then
      [xExtr, yErr1] = xMinimum(coef, xExtr, yErr1);
    else //  if (yErr1 < 0) then
      [xExtr, yErr1] = xMaximum(coef, xExtr, yErr1);
      end;
    
    jj = jj + 1;
    xx(jj) = xExtr;
    disp(string(xx(jj)) + "  " + string(yErr1));
    end; // while
  
  
  // Détermination du polynôme d'interpolation
  // Système d'équations à résoudre.
  ay = fct(xx);
  A = ones(ay);
  for ii=1:nnn;  A = [A, xx.^ii]; end;
  for ii=1:nnn+2;  A(ii, nnn+2) = (-1)^ii; end;
  
  // Coefficient du polynome passant par les nnn+1 points défini par ax, ay
  coef = A\ay;
  disp("erreur max env. = " + string(coef($)));
  err_max = coef($);
  coef = coef(1:$-1);
end; // for iter=...

scf(4);
clf()
plot(xxx, poly_f(coef, xxx) - fct(xxx));
title("Après une itération de Remez. Erreur de l''interpolation polynomiale de degré " + ..
      string(nnn));

// Comparaison entre l'erreur et le polynôme de Chebychev de degré n+1
plot(xxx, -err_max * horner(chepol(nnn+1, 'x'), 2*xxx-1), "g");

// Petite digression, sans grand intérêt
// Approximation meilleure sans trop de calculs
scf(5);
clf()
plot(xxx, poly_f(coef, xxx) + err_max * horner(chepol(nnn+1, 'x'), 2*xxx-1) - fct(xxx));
title("Amélioration gratuite par un polynôme de degré " + ..
      string(nnn+1));








