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

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

function ex21()
//=============
// Exercice 2.1.  Évolution d'une vitesse d'un corps qui chute avec frottement.
// Équation différentielle d'évolution de la vitesse.:
// V'(t) = 9.81 - (0.5/3.0) * V(t) , V(0) = 0 [m/s]
   function y_prime=f(t, y)
   //======================
   // t représente le temps  (non utilisé ici)
   // y  représente la vitesse
   // y_prime  représente la dérivée de la vitesse par rapport au temps
   g = 9.81; // [m/s^2]
   m = 3.0;  // masse du corps en [kg]
   kfrot = 1.5;   // coefficient de frottement [kg / s]
   y_prime = g - (kfrot/m)*y;  // car  vitesse ' = g - (kfrot/m) * vitesse
   endfunction
   
// données :
t0    =  0;   // temps initial
tFin  = 10;   // temps final
y0    = 0;  // vitesse initiale = 0 [m/s]
N  = 101;      // nombre de points d'évaluation
at = linspace(t0, tFin, N);  // vecteur LIGNE

// 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(1);  clf();
plot(at, ay, 'b');  // 'b' = blue, pour les positions

gcf().figure_name="Évolution de la vitesse en fonction du temps" // Titre de la fenêtre
title("Évolution de la vitesse en fonction du temps", 'fontsize', 3);
xlabel("Temps [s]", 'fontsize', 3);
ylabel("Vitesse [m/s]", 'fontsize', 3);
endfunction
//ex21();

function ex31()
//=============
// Exercice 3.  Évolution d'une vitesse d'un corps qui chute avec frottement.
// Avec l'évolution de la position en fonction du temps.
// Équation différentielle d'évolution de la vitesse.:
// X'(t) = V(t) , X(0) = 0 [m]
// V'(t) = 9.81 - (0.5/3.0) * V(t) , V(0) = 0 [m/s]

   function y_prime=f(t, y)
   //======================
  // t représente le temps  (non utilisé ici)
  // y(1)  représente la position
  // y(2)  représente la vitesse
  // y_prime(1)  représente la dérivée de la position par rapport au temps
  // y_prime(2)  représente la dérivée de la vitesse par rapport au temps
   // y et y_prime sont des vecteurs colonne
   g = 9.81; // [m/s^2]
   m = 3.0;  // masse du corps en [kg]
   kfrot = 1.5;   // coefficient de frottement [kg / s]
   V = y(2)
   y_prime = [ V; ..   // car  position' = vitesse
               (g - (kfrot/m)*V)];  // car  vitesse ' = g - (kfrot/m) * vitesse
   // Autre écriture possible
   //y_prime(1) = y(2);     // car  position' = vitesse
   //y_prime(2) = g - (kfrot/m)*y(2);  // car  vitesse ' = g - (kfrot/m) * vitesse
   endfunction
   
// données :
t0    =  0;   // temps initial
tFin  = 10;   // temps final
V0 = 0;       // vitesse initiale = 0 [m/s]
X0 = 0;       // position = 0 [m]
y0    = [X0; V0];  // vecteur COLONNE
N  = 101;      // nombre de points d'évaluation
at = linspace(t0, tFin, N);  // vecteur LIGNE

// 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(1);  clf();
plot(at, ay(2,:), 'b');  // 'b' = blue, pour les vitesses

gcf().figure_name="Évolution de la vitesse en fonction du temps" // Titre de la fenêtre
title("Évolution de la vitesse en fonction du temps", 'fontsize', 3);
xlabel("Temps [s]", 'fontsize', 3);
ylabel("Vitesse [m/s]", 'fontsize', 3);

scf(2);  clf();
plot(at, ay(1,:), 'g');  // 'g' = red, pour les positions

gcf().figure_name="Évolution de la position en fonction du temps" // Titre de la fenêtre
title("Évolution de la position en fonction du temps", 'fontsize', 3);
xlabel("Temps [s]", 'fontsize', 3);
ylabel("Position [m]", 'fontsize', 3);
endfunction
//ex31();

function ex42()
//=============
// Exercice 4
// problème 2 de la série 5.  Évolution du niveau d'eau d'un réservoir
// Équation différentielle d'évolution de la vitesse.:
// h'(t) = 2 - h/t) / 6 ,  h(0) = 0 [m]
// Il est clair que le débit D(t) = h'(t) sera nul si h(t) = 12 dm et
// même négatif au dela de 12 dm.  
// Donc la hauteur maximal est de 12 dm.
   function y_prime=f(t, y)
   //======================
  // t représente le temps  (non utilisé ici)
   // y  représente la hauteur  h
   // y_prime  représente  h',  la dérivée de la hauteur par rapport au temps

   y_prime = 2 - y / 6;  // car  h'(t) = 2 - h/t) / 6
   endfunction
   
// données :
t0    =  0;   // temps initial
tFin  = 30;   // temps final
y0    = 0;  // hauteur initiale = 0 [m/s]
N  = 101;      // nombre de points d'évaluation
at = linspace(t0, tFin, N);  // vecteur LIGNE

// 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(3);  clf();
plot(at, ay, 'b');  // 'b' = blue, pour la hauteur

gcf().figure_name="Évolution de la hauteur en fonction du temps" // Titre de la fenêtre
title("Évolution de la hauteur en fonction du temps", 'fontsize', 3);
xlabel("Temps [s]", 'fontsize', 3);
ylabel("hauteur [m]", 'fontsize', 3);
endfunction
//ex42();

function ex43()
//=============
// Exercice 3.  Évolution d'un oscillateur amorti
// 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)

  // Définition de l'équation différentielle :  [y(1) ; y(2)]' = [y(2) ; -k/m * y(1) -c/m * y(2)]
  function y_prime=f(t, y)
  //======================
  // t représente le temps  (non utilisé ici)
  // y(1)  représente la position
  // y(2)  représente la vitesse
  // y_prime(1)  représente la dérivée de la position par rapport au temps
  // y_prime(2)  représente la dérivée de la vitesse par rapport au temps
  // y et y_prime sont des vecteurs colonne
  m = 0.80;  // masse du corps
  k = 2.0;
  c = 0.30;
  X = y(1);
  V = y(2);
 
  y_prime = [ V, ..  // car  position' = vitesse
              (-k/m * X  - c/m * V)];  // car  vitesse' = -k/m * position - c/m * vitesse
  // Autre écriture possible
  //y_prime(1) = y(2);     // car  position' = vitesse
  //y_prime(2) = -k/m * y(1) -c/m * y(2);  // car  vitesse' = -k/m * position - c/m * vitesse
  endfunction

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

// 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,:), 'b');  // 'b' = blue, pour les positions
//plot(at, ay(2,:), 'r');  // 'r' = red , pour les vitesses
endfunction
//ex43();

