// Série 05,  résolution numérique de systèmes d'équations différentielles.

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

// Exercice 3 de la série 4.  Évolution d'un oscillateur avec frottement
// Décrit par une équation différentielle du second ordre :
// x(t)'' = -k/m * x(t) -c/m * v(t))
// ou par un système de deux équations différentielles du première ordre :
// x'(t) = v(t)
// v'(t) = -k/m * x(t) -c/m * v(t)
function ex3a()
//=============
// données :
m = 0.80;  // masse du corps
k = 2.0;   // constante de rappel
c = 0.3;   // Coef. de frottement proportionnelle à la vitesse
N  = 800;    // nombre de points d'évaluation
t0    = 0;    // temps initial
tFin = 20.0; // temps en secondes
dt = tFin / N;  // pas de la méthode numérique de résolution d'équ. diff.
at(1) = 0.0; // temps initial
ax(1) = 0.1; // Position initiale
av(1) = 0.0; // Vitesse initiale

// Résolution numérique
j = 0;  // indice pour mémoriser l'évolution
while (j < N)
   j = j + 1;
   at(j+1) = j*dt;
   a_de_t = -(k * ax(j) + c * av(j)) / m;  // accélération au temps (t-dt)
   ax(j+1) = ax(j) + dt * av(j);  // évolution de la position
   av(j+1) = av(j) + dt * a_de_t; // évolution de la vitesse
end;

scf(1);
clf();
plot(at, ax, 'k'); 
endfunction

// Exercice 3b.  Évolution d'un oscillateur avec frottement
// Décrit par une équation différentielle du second ordre :
// x(t)'' = -k * x(t) -c * v(t))
// ou par un système de deux équations différentielles du première ordre :
// x'(t) = v(t)
// v'(t) = -k/m * x(t) -c/m * v(t)
function ex3b()
//============
// Utilise la méthode de Runge
// données :
m = 0.80;  // masse du corps
k = 2.0;   // constante de rappel
c = 0.3;  // Coef. de frottement proportionnelle à la vitesse
N  = 80;      // nombre de points d'évaluation
t0    = 0;    // temps initial
tFin = 20.0; // temps en secondes
dt = tFin / N;  // pas de la méthode numérique de résolution d'équ. diff.
at(1) = 0.0; // temps initial
ax(1) = 0.1; // Position initiale
av(1) = 0.0; // Vitesse initiale

// Résolution numérique
j = 0;  // indice pour mémoriser l'évolution
while (j < N)
   j = j + 1;
   at(j+1) = j*dt;
   a_de_t = -(k * ax(j) + c * av(j)) / m;  // accélération au temps (t-dt)
   x1s2 = ax(j) + dt/2 * av(j);  // position en (t - dt/2)
   v1s2 = av(j) + dt/2 * a_de_t; // vitesse  en (t - dt/2)
   a_de_t_plus_dts2 = -(k*x1s2 + c*v1s2 ) / m;
   ax(j+1) = ax(j) + dt * v1s2;
   av(j+1) = av(j) + dt * a_de_t_plus_dts2;
end;

// graphique de la solution numérique
//scf(2);  clf();
plot(at, ax(:), 'b');  // axv(1, :) est un vecteur ligne.
endfunction

// Exercice 3c.  Évolution d'un oscillateur avec frottement
// Décrit par une équation différentielle du second ordre :
// x(t)'' = -k * x(t) -c * v(t))
// ou par un système de deux équations différentielles du première ordre :
// x'(t) = v(t)
// v'(t) = -k/m * x(t) -c/m * v(t)
function ex3c()
//============
// Utilise la méthode de Runge
// Fonction écrite sous forme vectorielle
   function z=f(y)
   //=============
   // y(1)  représente la position
   // y(2)  représente la vitesse
   // z(1)  représente la dérivée de la position par rapport au temps
   // z(2)  représente la dérivée de la vitesse par rapport au temps
   // y et z sont des vecteurs colonne
   m = 0.80;  // masse du corps
   k = 2.0;   // constante de rappel
   c = 0.3;  // Coef. de frottement proportionnelle à la vitesse
   z(1) = y(2);     // car  position ' = vitesse
   z(2) = -k/m * y(1) -c/m * y(2);  // car  vitesse ' = -k * position - c * vitesse
   endfunction
   
// données :
N  = 50;      // nombre de points d'évaluation
t0    = 0;    // temps initial
tFin = 20.0; // temps en secondes
dt = tFin / N;  // pas de la méthode numérique de résolution d'équ. diff.
at(1) = 0.0; // temps initial
axv(:, 1) = [0.1, 0.0]; // Position initiale et Vitesse initiale
   
// Résolution numérique
j = 0;  // indice pour mémoriser l'évolution
while (j < N)
   j = j + 1;
   at(j+1) = j*dt;
   k1 = f(axv(:, j));
   k2 = f(axv(:, j) + k1 * dt/2);
   k3 = f(axv(:, j) + k2 * dt/2);
   k4 = f(axv(:, j) + k3 * dt);   
   axv(:, j+1) = axv(:, j) + (k1 + 2*k2 + 2*k3 + k4)*dt/6;
end;

// graphique de la solution numérique
//scf(3);  clf();
plot(at', axv(1, :), 'r');
endfunction

global cpt;  // compte le nombre d'appels à la fonction  f

// Résolution du même système, avec la méthode intégrée dans Scilab
function ex3d()
//=============
  // Définition de l'équation différentielle :  [y(1) ; y(2)]' = [y(2) ; -k/m * y(1) -c/m * y(2)]
  function z=f(t, y)
  //================
  // t représente le temps  (non utilisé ici)
  // y(1)  représente la position
  // y(2)  représente la vitesse
  // z(1)  représente la dérivée de la position par rapport au temps
  // z(2)  représente la dérivée de la vitesse par rapport au temps
  m = 0.80;  // masse du corps
  k = 2.0;
  c = 0.30;
  z(1) = y(2);     // car  position' = vitesse
  z(2) = -k/m * y(1) -c/m * y(2);  // car  vitesse' = -k/m * position - c/m * vitesse
  global cpt;
  cpt = cpt+1;
  endfunction

// données :
t0    =  0;   // temps initial
tFin  = 20;   // temps final
y0    = [0.1; 0];  // position = 0.1 [m] ; vitesse initiale = 0 [m/s]], vecteur COLONNE
N  = 101;      // nombre de points d'évaluation
at = linspace(t0, tFin, N);  // vecteur LIGNE

global cpt;
cpt = 0; // compte le nombre d'appels à la fonction  f   
// résolution numérique en appelant la fonction  ode  de SciLab
ay = ode( y0,..     // condition initiale
          t0,..     // t_initial,
          at,..     // série des t(j) pour lesquels  y(j)  sera calculé 
          1e-3,..   // tolérance d'erreur relative
          1e-4,..   // tolérance d'erreur absolue
          f..       // fonction définissant l'équation différentielle
          );
// L'algorithme essaye d'avoir une erreur inférieure au
// maximum entre la tolérance d'erreur absolue et la tolérance d'erreur relative.

// graphique de la solution numérique
//scf(4);  clf();
plot(at, ay(1,:), 'g');  // 'b' = blue, pour les positions
//plot(at, ay(2,:), 'r');  // 'r' = red , pour les vitesses
disp('Nombre d''appels à la fonction  f : ' + string(cpt));
endfunction

// Exécute les 4 fonctions
ex3a();  ex3b();  ex3c();  ex3d();
