// Approximation par les polynôme de Cebyshev (Tchebychev)  Tk(x)
// A comparer avec le programme : s04c_approximation_polynomiale.sce
// Il y a peut-être une amélioration de précision de 1 bit,
// lorsqu'on atteint les limites de précision de la machine.
// Sinon, on ne remarque pas de différence avec le cas
// d'utilisation de polynôme écrit de manière standard.
// C.f. le livre : 
// "Numerical recipes in C, the Art of Scientific Computing"
// Edition Cambridge
// "Chebyshev approximation"  chapitre 5

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+1) / sqrt(2))/2; // à partir de nnn=18, on atteind la précision de l'ordinateur
y = sin((x+1)*%pi/4); // à partir de nnn=14, on atteind la précision de l'ordinateur
//y = sin(sqrt(x+1)*%pi/4) ./ sqrt(x+1);
//y = exp(x+1); // excellante approximation de degrés 4, 6, 8.  nnn=10 atteint les limites de précision
//y = cos(x*%pi/4); // pour nnn > 10, les limites de précision sont atteintes
//y = sqrt(2.5 + 1.5*x);
//y = 1 ./ (1 + 25*(x+1).^2); // Dans aucun cas l'approximation est excellante !
endfunction

// Degré du polynôme d'interpolation = nombre de points - 1 = nnn
nnn = 8;

// Points à interpoler :
aaa=-1;  // intervalle d'interpolation
bbb= 1;
xxx = linspace(aaa+1e-16, bbb, 20*nnn+1);  // pour les plot
kkk = (0:nnn)';
ax = cos((2*kkk+1)*%pi / (2*nnn+2));
//ax = linspace(aaa, bbb, nnn+1)';  // moins bonne approximation
ay = fct(ax);

// Evaluation du polynôme, où  coef  sont les coefficients des pol. de Cebyshev
function y=poly_f(coef, x)
//========================
// Polynome de degré  nnn, de coefficients coef
djp2 = zeros(x);
djp1 = zeros(x);
for j=length(coef)-1:-1:1 
  dj = 2 * x .* djp1 - djp2 + coef(j+1) * ones(x);
  djp2 = djp1;
  djp1 = dj;
  end;
y = x .* djp1 - djp2 + coef(1) * ones(x);  // ici  j = 0;
endfunction

scf(0);
clf(0);
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];
title("Courbe à approximer.");

// Détermination du polynôme d'interpolation passant par les points de ax, ay
// Système d'équations à résoudre.
A = [];
for ii=0:nnn;  A = [A, horner(chepol(ii, 'x'), ax)]; end;

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

// autre manière de calculer les coefficients,
// lorsque :  ax = cos((2*kkk+1)*%pi / (2*nnn+2));  Pas valide sinon.
coef3 = zeros(ax);
for j=0:nnn
  coef3(j+1) = (2/(nnn+1)) * sum(fct(ax) .* cos( j*%pi*(2*kkk+1) / (2*nnn+2)));
  end;
coef3(1) = 0.5*coef3(1);
// coef3 == coef  à la précision de la machine près
// abs(coef - coef3 ) <= 2*%eps   donne une série de TRUE

scf(1);
clf();
plot(xxx, poly_f(coef, xxx) - fct(xxx));
title("Différence entre le polynôme d''interpolation de degré " +..
      string(length(ax)-1) + " et la fonction d''origine");

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

// MIEUX, légèrement !
// Détermination du polynôme d'interpolation ayant les extremums 
// des différences avec la fonction d'origine en  axx  qui alternent
// Système d'équations à résoudre.
kkkk = (0:nnn+1)';
axx = cos((kkkk*%pi / (nnn+1))) + 1e-16;  // +1e-16 pour éviter les div. par 0
ayy = fct(axx);
AA = [];
for ii=0:nnn;  AA = [AA, horner(chepol(ii, 'x'), axx)]; end;
for ii=1:nnn+2;  AA(ii, nnn+2) = (-1)^ii; end;

// Coefficient du polynome passant par les nnn+1 points défini par ax, ay
coef2 = AA\ayy;   // coef2($) env. = l'erreur max d'approximation

scf(2);
clf();
plot(xxx, poly_f(coef2(1:$-1), xxx) - fct(xxx));
title("Différence entre le polynôme d''interpolation de degré " +..
      string(nnn) + " et la fonction d''origine");
