// Série 10, divers problèmes pour s'exercer au questions de maturité.

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


function coef = FourCoef(af)
//=========================
// Retourne dans rep(1, 2:M+1) les coefficients  Ak  de Fourier
// Retourne dans rep(2, 2:M+1) les coefficients  Bk  de Fourier
// rep(1, 1) = A0
// rep(2, 1) =  0
N = size(af, 2);
M = N / 2 - 1;  //disp("M=" + string(M));

A0 = sum(af) / N;

at = 2*%pi/N*(0:N-1)'; // Pour avoir un vecteur colonne
aA = zeros(1,M);   aB = zeros(1,M);  // Pour imposer d'avoir des vecteurs lignes.
for k=1:M do
  aA(k) = af * cos(k*at) * 2 / N
  aB(k) = af * sin(k*at) * 2 / N
end

coef = [A0, aA ; 0, aB];
endfunction

function coef = FourCoef2(af)
//==========================
// Version moins efficace, mais plus intuitive à programmer
// Retourne dans rep(1, 2:M+1) les coefficients  Ak  de Fourier
// Retourne dans rep(2, 2:M+1) les coefficients  Bk  de Fourier
// rep(1, 1) = A0
// rep(2, 1) =  0
N = size(af, 2);
M = N / 2 - 1;  //disp("M=" + string(M));

A0 = 0;
for j=1:N do
    A0 = A0 + af(j)
end
A0 = A0 / N;

at = 2*%pi/N*(0:N-1)'; // Pour avoir un vecteur colonne
aA = zeros(1,M);   aB = zeros(1,M);  // Pour imposer d'avoir des vecteurs lignes.

// Boucle pour définir chaque coefficient aA(k)  et  aB(k)
for k=1:M do
  aA(k) = 0; // Initialisation à 0
  aB(k) = 0;
  // Boucle pour effectuer la somme désirée.
  for j=1:N do
    aA(k) = aA(k) + af(j) * cos(k*at(j)) * 2 / N
    aB(k) = aB(k) + af(j) * sin(k*at(j)) * 2 / N
  end
end

coef = [A0, aA ; 0, aB];
endfunction

function ex11()
//=============
// Recherche des fréquences et amplitudes qui caractérisent une fonction périodique.
// Signal a  w1/(2 pi)  et  w2/(2 pi)  Hz
w1 =   5;  A1 = 1.3;
w2 =   15;  B2 = 0.6;
N = 2^7;
at = 2*%pi/N*(0:N-1)
fonc  = A1*cos(w1*at) + B2*sin(w2*at);

scf(11);
clf();
plot(at, fonc, 'k')

Coefs = FourCoef2(fonc)

win2 = scf(12);
clf();
plot(Coefs(1, 2:23), 'g');
plot(Coefs(2, 2:23), 'r');

axes2 = win2.children;
title( axes2, "Coefficients A(nu) et B(nu)", 'fontsize', 4);
xlabel(axes2, "Pulsation  ''w''  en [rad / s]", 'fontsize', 3);
ylabel(axes2, "Amplitudes", 'fontsize', 3);
endfunction
//ex11();

function ex12()
//=============
// Recherche des fréquences et amplitudes qui caractérisent une fonction périodique.

// Pour indiquer le dossier de travail
realpath = cd("/home/bg/Bg/Cours/AppMath_4/Exercices");

disp(listvarinfile("s10d_ex1a_data.bin"));

// Pour charger des données au format de Scilab
// Les données sont sauvegardées dans la variable "fonc"
// Le nom de la variable "fonc" est obligatoire.
load('s10d_ex1a_data.bin', 'fonc');

N = size(fonc, 2);
at = 2*%pi/N*(0:N-1)

scf(13);
clf();
plot(at, fonc, 'k')

Coefs = FourCoef(fonc)

win2 = scf(14);
clf();
plot(Coefs(1, 2:25), 'g');
plot(Coefs(2, 2:25), 'r');

axes2 = win2.children;
title( axes2, "Coefficients A(nu) et B(nu)", 'fontsize', 4);
xlabel(axes2, "Pulsation  ''w''  en [rad / s]", 'fontsize', 3);
ylabel(axes2, "Amplitudes", 'fontsize', 3);
endfunction
//ex12();


// Pour la préparation de l'exercice.
function ex10()
//=============
// Génération d'un fichier contenant un signal périodique de période 2 pi
// L'objectif est de déterminer les fréquences et amplitudes ce ce signal.
w1 =   7;   A1 = 1.7;
w2 =  12;   B2 = 0.8;
w3 =  19;   A3 = 1.2;

N = 2^7;
at = 2*%pi/N*(0:N-1)
fonc  = A1*cos(w1*at) + B2*sin(w2*at) + A3*cos(w3*at);

// Pour indiquer le dossier de travail
realpath = cd("/home/bg/Bg/Cours/AppMath_4/Exercices");

// Pour sauvegarder des données au format de Scilab
// Ce sont les données du vecteur "fonc" qui sont sauvegardés.
save("s10d_ex1a_data.bin", "fonc");

// Ce qui suit est juste pour des vérifications
disp(listvarinfile("s10d_ex1a_data.bin"));

scf(1);
clf();
plot(at, fonc, 'k')

Coefs = FourCoef(fonc)

win2 = scf(2);
clf();
plot(Coefs(1, 2:23), 'g');
plot(Coefs(2, 2:23), 'r');

axes2 = win2.children;
title( axes2, "Coefficients A(nu) et B(nu)", 'fontsize', 4);
xlabel(axes2, "Pulsation  ''w''  en [rad / s]", 'fontsize', 3);
ylabel(axes2, "Amplitudes", 'fontsize', 3);
endfunction
//ex10();

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

function ex21()
//=============
// Détermination de la capacité d'un condensateur,
// basé sur des mesures de décharge du condensateur à travers une résistances fonction du temps.
// Basée sur le travail de Arthur Freeman, Date: 07/02/2019
// Avec les données expérimentales de Arthur Freeman

R = 1000; // [Ohms], la valeur de la résistance

// Vecteur des mesures de temps
mesuresTemps = [    0,  5.70,  11.48,  18.31,  25.31,  32.65,  41.02,  49.46,  59.02,  69.51, ..
                81.17, 93.74, 108.22, 124.31, 143.71, 166.70, 195.34, 233.25, 288.88, 394.31]'
// Vecteur des tensions de mesures
tensions = linspace(10, 0.5, 20)';

scf(11); clf()//On initialise la fenêtre graphique 11.
plot(mesuresTemps, tensions, 'x', 'thickness', 1),xgrid 

title('Tension en fonction du temps.','fontsize',5, 'fontname', 'helvetica')
ylabel('Uc(t) [V]','fontsize',4) //Esthétique.
xlabel('t [s]','fontsize',4)

//===========================================================
// Logarithme naturel des tensions
log_tensions = log(tensions)

scf(12); clf()//On initialise la fenêtre graphique
plot(mesuresTemps, log_tensions, 'x', 'thickness', 1),xgrid 

// Matrice de détermination des deux inconnues :
// Log de la tension initiale
// Produit -1/(R*C)
A  = ones(log_tensions);
A = [A, mesuresTemps]; // matrice du système d'équations

coef = A \ log_tensions;  // coefficients de la fonction cherchée
disp("Tension initiale = " + string(exp(coef(1))))

C_fit_lineaire = -1 / (R*coef(2));
disp("Capacité = " + string(C_fit_lineaire))

// Pour avoir la courbe déterminée par régression linéaire
log_analytique = coef(1) + coef(2) * mesuresTemps;
plot( mesuresTemps, log_analytique, 'r');

title('Log(Tension) en fonction du temps.','fontsize',4, 'fontname', 'helvetica')
ylabel('log(Uc(t))','fontsize',4) //Esthétique.
xlabel('t [s]','fontsize',4)

// Dessin de la courbe obtenue par régression linéaire, dans le graphique de la décharge
scf(11);
analytique = exp(coef(1)) * exp(coef(2) * mesuresTemps);
plot( mesuresTemps, analytique, 'r');

//===========================================================
// On se limite au 15 premières mesures, car les suivantes sont moins fiables,
// car le condensateur est presque déchargé.
// Les résultats obtenus sont bien meilleurs en ayant éliminé les 5 dernières mesures !
tensions = tensions(1:15);
mesuresTemps = mesuresTemps(1:15);

// Logarithme naturel des tensions
log_tensions = log(tensions);

// Matrice de détermination des deux inconnues :
// Log de la tension initiale
// Produit -1/(R*C)
A  = ones(log_tensions);
A = [A, mesuresTemps]; // matrice du système d'équations

coef = A \ log_tensions;  // coefficients de la fonction cherchée
disp("Tension initiale = " + string(exp(coef(1))))

C_fit_lineaire = -1 / (R*coef(2));
disp("Capacité = " + string(C_fit_lineaire))

// Pour avoir la courbe déterminée par régression linéaire
log_analytique = coef(1) + coef(2) * mesuresTemps;
scf(12)
plot( mesuresTemps, log_analytique, 'k');

// Dessin de la courbe obtenue par régression linéaire, dans le graphique de la décharge
scf(11);
analytique = exp(coef(1)) * exp(coef(2) * mesuresTemps);
plot( mesuresTemps, analytique, 'k');
endfunction
//ex21();

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

function ex31()
//=============
// Approximation de la fonction racine carrée sur l'intervalle [1 ; 4] par des
// polynomes de Tchebychev.   Tk(x) = cos(k*arccos(x)), k=0,1,2,...
// On veut que  f(x_j) = a1 * T_0(x_j) + a2 * T_1(x_j) + a3 * T_2(x_j) + ...  an * T_(n-1)(x_j).  n = 7
// f(x) = racine(1.5*x + 2.5)  x  dans  [-1 ; 1]
// C'est pas évident, mais les polynômes de Tchebychev sont bien des polynômes.
// Ils satisfont aussi la propriété :  T_k(x) = 2*x*T_k-1(x) - T_k-2(x),  k >= 2.
//
// Pour obtenir une meilleure approximation, il faut augmenter la valeur de  n.

// Points par lesquels la courbe doit passer au mieux. On en prend une quinzaine.
ax = linspace(-1, 1, 11)';
ay = sqrt(1.5*ax + 2.5);  // Rapporte la racine carrée dans l'intervalle [-1 ; 1]

// On place les points dans un graphique
scf(31);  clf();
plot( ax, ay, 'LineStyle', 'none', "Marker", 'o');

// Détermination de la matrice qui permet de trouver les coefficients  ak  de l'approximation
n = 5;
A  = ones(ax);
for kk=2:n do
  A = [A, cos((kk-1)*acos(ax))]; // matrice du système d'équations
  end;
ak = A\ay;  // coefficients de la fonction cherchée

// Graphique de l'approximation et de la fonction d'origine à approximer
xxx = linspace(-1, 1, 101)';
Fxxx = zeros(xxx);
for kk=1:length(ak) do
   Fxxx = Fxxx + ak(kk)*cos((kk-1)*acos(xxx));
   end;

// On trace la courbe qui passe au mieux par les points
plot(xxx, Fxxx); 

// différence entre les points et la courbe
scf(32);  clf();
aFx = zeros(ax);
for kk=1:length(ak) do
   aFx = aFx + ak(kk)*cos((kk-1)*acos(ax));
   end;
plot(ax, aFx - ay, 'thickness', 2); // On trace la courbe

// -------------------------------------------------------------------------
// La même approximation avec un polynôme au lieu des fonctions cos(k*acos(x))
// On veut donc que :  f(x_j) = a1 + a2 * x_j + a3 * (x_j)^2 + ...  an * (x_j)^(n-1).  n = 7
// f(x) = racine(1.5*x + 2.5)  x  dans  [-1 ; 1]
// On remarque que l'on obtient exactement les mêmes approximations !
A  = ones(ax);
for kk=2:n do
  A = [A, ax.^(kk-1)]; // matrice du système d'équations
  end;
ak = A\ay;  // coefficients de la fonction cherchée

scf(31);
plot( ax, ay, 'LineStyle', 'none', "Marker", 'o'); // On place les points

// Graphique de l'approximation et de la fonction d'origine à approximer
xxx = linspace(-1, 1, 101)';
Fxxx = zeros(xxx);
for kk=1:length(ak) do
   Fxxx = Fxxx + ak(kk)*xxx.^(kk-1);
   end;

plot(xxx, Fxxx, 'LineStyle', '-.', 'Color', 'red'); // On trace la courbe

scf(32); 
aFx = zeros(ax);
for kk=1:length(ak) do
   aFx = aFx + ak(kk)*ax.^(kk-1);
   end;
plot(ax, aFx - ay, 'LineStyle', '-.', 'Color', 'red', 'thickness', 2); // On trace la courbe
endfunction
//ex31();






