// Exercices de la série 4 du cours d'OS AM : matrices, interpolation et approximation

funcprot(0); // pour ne pas afficher le message : "redefining function".


function ex01()
//=============  
// Détermination de la droite passant au mieux par deux points données
ax = [-2; 3];  // Vecteur colonne
ay = [1 ; 5];

// Matrice du système d'équations à résoudre.
A = [1, ax(1) ;..
     1, ax(2)];

// Coefficient  b  et  a  de la droite d'équation  y = a*x + b passant par les deux points définis par ax, ay
coef = A\ay;
disp(coef, "coef =");

win = scf(1);  // Référence à la fenêtre contenant le graphique.
clf(1); // efface le contenu de la fenêtre.

// Référence au graphique
axes = win.children;  // la courbe et les axes contenu dans la fenêtre graphique

// Affiche les 2 points par lesquels passent le polynôme.
plot(ax, ay, 'LineStyle', 'none', 'Marker', 'o');

// Points pour tracer la courbe de la droite
xxx = linspace(min(ax)-1, max(ax)+1, 2);

// Polynôme de degré 1 passant par les 2 points définis par les vecteurs  ax  et  ay.
yyy =  coef(1) + coef(2) .* xxx;

// tracé du graphique.
plot(xxx, yyy);

// Définition de l'échelle des axes
axes.zoom_box = [min(ax)-0.2, min(ay)-0.2, max(ax)+0.2, max(ay)+0.2];
title("Droite passant par les deux points donnés", 'fontsize', 3);
xlabel(axes, "axe des x", 'fontsize', 3);
ylabel(axes, "axe des y", 'fontsize', 3);
endfunction
//ex01();


function ex02()
//=============  
// Détermination de la droite passant au mieux par deux points données
ax = [-2; 3; 8; 13];  // Vecteur colonne
ay = [1 ; 5; 8; 13];

// Matrice du système d'équations à résoudre.
A = [1, ax(1) ;..
     1, ax(2) ;..
     1, ax(3) ;..
     1, ax(4) ;..
     ];

// Coefficient  b  et  a  de la droite d'équation  y = a*x + b passant 
// au mieux par les points définis par ax, ay
coef = A\ay;
disp(coef, "coef =");

win = scf(1);  // Référence à la fenêtre contenant le graphique.
clf(1); // efface le contenu de la fenêtre.

// Référence au graphique
axes = win.children;  // la courbe et les axes contenu dans la fenêtre graphique

// Affiche les points par lesquels passent le polynôme.
plot(ax, ay, 'LineStyle', 'none', 'Marker', 'o');

// Points pour tracer la courbe de la droite
xxx = linspace(min(ax)-1, max(ax)+1, 2);

// Polynôme de degré 1 passant au mieux par les points définis par les vecteurs  ax  et  ay.
yyy =  coef(1) + coef(2) .* xxx;

// tracé du graphique.
plot(xxx, yyy);

// Définition de l'échelle des axes
axes.zoom_box = [min(ax)-0.5, min(ay)-0.5, max(ax)+0.5, max(ay)+0.5];
title("Droite passant au mieux par les points donnés", 'fontsize', 3);
xlabel(axes, "axe des x", 'fontsize', 3);
ylabel(axes, "axe des y", 'fontsize', 3);
endfunction
//ex02();


function ex03()
//=============  
// Détermination du polynôme d'interpolation passant par les points donnés
ax = [-2; -1; 0; 1; 2];  // Vecteur colonne
ay = sin(ax);

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

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

win = scf(1);  // Référence à la fenêtre contenant le graphique.
clf(1); // efface le contenu de la fenêtre.

// Référence au graphique
axes = win.children;  // la courbe et les axes contenu dans la fenêtre graphique

// Affiche les points par lesquels passent le polynôme.
plot(ax, ay, 'LineStyle', 'none', 'Marker', 'o');

// Points pour tracer la courbe du polynôme.
xxx = linspace(min(ax), max(ax), 101);

// Polynôme de degré 4 passant par les 5 points définis par les vecteurs  ax  et  ay.
yyy =  coef(1) + coef(2) .* xxx + coef(3) .* xxx.^2 + coef(4) .* xxx.^3 + coef(5) .* xxx.^4;

// tracé du graphique.
plot(xxx, yyy);

// Définition de l'échelle des axes
axes.zoom_box = [min(ax)-0.1, -1.1, max(ax)+0.1, 1.1];
title("Polynôme d''interpolation de degré " + string(length(ax)-1));
xlabel(axes, "axe des x", 'fontsize', 3);
ylabel(axes, "axe des y", 'fontsize', 3);

endfunction
//ex03();


function ex04()
//=============  
// Détermination du polynôme d'interpolation passant par les points donnés
ax = linspace(-2, 2, 21)';  // Vecteur colonne
ay = sin(ax);

// Degré du polynôme passant au mieux par les points donnés.
nDeg = 3;

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

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

win = scf(1);  // Référence à la fenêtre contenant le graphique.
clf(1); // efface le contenu de la fenêtre.

// Référence au graphique
axes = win.children;  // la courbe et les axes contenu dans la fenêtre graphique

// Affiche les points par lesquels passent le polynôme.
plot(ax, ay, 'LineStyle', 'none', 'Marker', 'o');

// Points pour tracer la courbe du polynôme.
xxx = linspace(min(ax), max(ax), 101);

// Polynôme de degré  nDeg  passant au mieux par les points définis par les vecteurs  ax  et  ay.
yyy = coef(1) * ones(xxx);
for kk=1:nDeg
  yyy =  yyy + coef(kk+1) .* xxx.^kk;
  end

// tracé du graphique.
plot(xxx, yyy);

// Définition de l'échelle des axes
axes.zoom_box = [min(ax)-0.1, -1.1, max(ax)+0.1, 1.1];
title("Polynôme de degré " + string(nDeg) + " passant au mieux par les points donnés");
xlabel(axes, "axe des x", 'fontsize', 3);
ylabel(axes, "axe des y", 'fontsize', 3);

endfunction
//ex04();


function ex5()
//============
// Approximation de la fonction exp(-alpha*x^2) par une somme de cosinus

  function y=Fourier(coefs, x)
  //==========================
  // y = coef(1) + coef(2)*cos(x) + coef(3)*cos(2*x) + ... + coef(n)*cos((n-1)*x)  n=length(coefs)
  y = zeros(x);
  for ii=1:length(coefs) do
     y = y + coefs(ii)*cos((ii-1)*x);
     end;
  endfunction
  
Periode = 2*%pi;
alpha = 5;
ax=linspace(-Periode/2, Periode/2, 15)';
ay = exp(-alpha*ax.^2);
A = [ones(ax)];
for ii=2:7 do
  A = [A, cos((ii-1)*ax)]; // matrice du système d'équations
  end;
coefs = A\ay;  // coefficients de la fonction cherchée
//disp("Coefficients du terme constant au facteur de cos((n-1)*x).")
disp(coefs);

win0 = scf(1);
clf(1); 
plot( ax, ay, 'LineStyle', 'none', "Marker", 'o'); // On place les points
beta = 1;
ttt = linspace(beta*ax(1), beta*ax($), 101)';
plot(ttt, Fourier(coefs, ttt)); // On trace la courbe
plot(ttt, exp(-alpha*ttt.^2), "r"); // On trace la courbe d'origine en rouge

axes0 = win0.children; // Pas demandé
title("exp(-" + string(alpha) + "*x^2) et une approximation de Fourier.", 'fontsize',3);
xlabel(axes0, "x ", 'fontsize',3);
ylabel(axes0, "exp(-" + string(alpha) + "*x^2) et Fourier(x)", 'fontsize',3);

// différences
win0 = scf(11);
clf(); 
plot(ttt, Fourier(coefs, ttt) - exp(-alpha*ttt.^2)); // On trace la courbe

axes0 = win0.children; // Pas demandé
title("Différence entre le exp(-" + string(alpha) + "*x^2) et une approximation de Fourier.",..
        'fontsize',3);
xlabel(axes0, "x [rad]", 'fontsize',3);
ylabel(axes0, "Fourier(x) - exp(-" + string(alpha) + " * x^2)", 'fontsize',3);
endfunction
//ex5();


