// Série 06,  résolution numérique d'équations différentielles

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

function ex1()
//============
  // Définition de l'équation différentielle :  
  // Simulation d'un circuit R C, c.f. ex. 1 de la série 6
  // Décharge d'un condensateur
  function z=f(t, y)
  //================
  // t représente le temps
  // y(1)  représente la tension aux bornes du condensateur    U  en volts
  // z(1)  représente la dérivée de y(1)
  C = 1.0e-4;  // valeur de la capacité en Farads
  R = 15e3;   // valeur de la résistance, en Ohms 
  z(1) = -y(1)/(R*C);  // UC' = - UC / RC
  endfunction

  // données :
  t0    =  0;   // temps initial
  tFin  =  5;   // temps final
  y0    = 12;  // UC(0) = 12 [V]
  N  = 500;      // nombre de points d'évaluation
  at = linspace(t0, tFin, N+1);  // vecteur LIGNE

  // résolution numérique en appelant la fonction  ode  de SciLab
  ay = ode( y0,..     // condition initiale
            t0,..     // x_initial,
            at,..     // série des t(j) pour lesquels  y(j)  sera calculé 
            1e-4,..   // tolérance d'erreur relative
            1e-5,..   // 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 l'erreur absolue et l'erreur relative.

  // graphique de la solution numérique
  win0 = scf(1);
  clf();
  plot(at, ay(1,:), 'b');  // 'b' = blue, pour les tensions aux bornes du condensateur
  axes0 = win0.children;
  title("Tension [V] en bleue", 'fontsize', 4);
  xlabel(axes0, "Temps [s]", 'fontsize', 3, "color", [0, 0, 0]);
  ylabel(axes0, "Tension [V]", 'fontsize', 3, "color", [0, 0, 1]);
endfunction
// ex1
// On remarque que la charge du condensateur vaut U0 / e = 4,4145 [V] après
// un temps égale à R * C = 1,5 secondes. 
// C'est bien ce qui est prédit par la théorie.

function ex2()
//============
  // Définition de l'équation différentielle :  
  // Simulation d'un circuit R C, c.f. ex. 2 de la série 6
  // Charge d'un condensateur
  function z=f(t, y)
  //================
  // t représente le temps
  // y(1)  représente la tension aux bornes du condensateur    U  en volts
  // z(1)  représente la dérivée de y(1)
  C = 1.0e-4;  // valeur de la capacité en Farads
  R = 15e3;   // valeur de la résistance, en Ohms 
  U0 = 12.0;  // volts
  z(1) = (U0 - y(1)) / (R*C);  // UC' = (U0 - UC) / RC   UC(t) = U0 * (1 - exp(-t/RC))
  endfunction

  // données :
  t0    =  0;   // temps initial
  tFin  =  5;   // temps final
  y0    =  0;  // UC(0) = 0 [V]
  N  = 500;      // nombre de points d'évaluation
  at = linspace(t0, tFin, N+1);  // vecteur LIGNE

  // résolution numérique en appelant la fonction  ode  de SciLab
  ay = ode( y0,..     // condition initiale
            t0,..     // x_initial,
            at,..     // série des t(j) pour lesquels  y(j)  sera calculé 
            1e-4,..   // tolérance d'erreur relative
            1e-5,..   // 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 l'erreur absolue et l'erreur relative.

  // graphique de la solution numérique
  win0 = scf(1);
  clf();
  plot(at, ay(1,:), 'b');  // 'b' = blue, pour les tensions aux bornes du condensateur
  axes0 = win0.children;
  title("Tension [V] en bleue", 'fontsize', 4);
  xlabel(axes0, "Temps [s]", 'fontsize', 3, "color", [0, 0, 0]);
  ylabel(axes0, "Tension [V]", 'fontsize', 3, "color", [0, 0, 1]);
endfunction
// ex2
// On remarque que la charge du condensateur vaut U0 * (1 - 1/e) = 7,585 [V] après
// un temps égale à R * C = 1,5 secondes. 
// C'est bien ce qui est prédit par la théorie.

function ex3()
//============
  // Définition de l'équation différentielle :  
  // Simulation d'un circuit R C, c.f. ex. 2 de la série 6
  // Charge d'un condensateur
      
  function z=f(t,y)
  //===============
  // t = le temps
  // y(1) = tension aux bornes du condensateur au temps  t.
  // Lors de la décharge : y'(1) = -y(1) / (R*C)
  // Lors de la charge :   y'(1) = (U0 - y(1)) / (R*C)
  // On charge si y(1) < Umin
  // On décharge si y(1) > Umax
  global fCharge;  // pour dire d'utiliser la variable fCharge déclarée à l'ext. de la fonction.
  
  C = 1.0e-4;  // valeur de la capacité en Farads
  R = 15e3;   // valeur de la résistance, en Ohms 
  U0 = 12.0;  // volts
  Umin =  1;  // volts
  Umax = 11;  // volts
  if (y(1) < Umin) fCharge = %t; end;
  if (y(1) > Umax) fCharge = %f; end;
  if (fCharge) z(1) = (U0 - y(1)) / (R*C); // cas de charge
  else         z(1) = ( 0 - y(1)) / (R*C); // cas de décharge
  end;  
  endfunction

  // données :
  fCharge = %t;
  t0    =  0;   // temps initial
  tFin  = 20;   // temps final
  y0    =  0;  // UC(0) = 0 [V]
  N  = 500;      // nombre de points d'évaluation
  at = linspace(t0, tFin, N+1);  // vecteur LIGNE

  // résolution numérique en appelant la fonction  ode  de SciLab
  ay = ode( y0,..     // condition initiale
            t0,..     // x_initial,
            at,..     // série des t(j) pour lesquels  y(j)  sera calculé 
            1e-4,..   // tolérance d'erreur relative
            1e-5,..   // 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 l'erreur absolue et l'erreur relative.

  // graphique de la solution numérique
  win0 = scf(1);
  clf();
  plot(at, ay(1,:), 'b');  // 'b' = blue, pour les tensions aux bornes du condensateur
  axes0 = win0.children;
  title("Tension [V] en bleue", 'fontsize', 4);
  xlabel(axes0, "Temps [s]", 'fontsize', 3, "color", [0, 0, 0]);
  ylabel(axes0, "Tension [V]", 'fontsize', 3, "color", [0, 0, 1]);
endfunction
// ex3

function ex4()
//============
  // Définition de l'équation différentielle :  
  // Simulation d'un circuit R L C, c.f. ex. 1 de la série 7
  function z=f(x, y)
  //================
  // x représente le temps
  // y(1)  représente la tension aux bornes du condensateur    UC  en volts
  // y(2)  représente l'intensité du courant traversant la self  IL  en ampères
  // z(1)  représente la dérivée de y(1)
  // z(2)  représente la dérivée de y(2)
  C = 1.0e-6;  // valeur de la capacité en Farads
  L = 4;  // valeur de la self (inductance) en Henri
  R = 15e3;   // valeur de la résistance, en Ohms 
  z(1) = (1/C) * (y(2) - y(1)/R);  // UC' = (1/C) * (IL - UC/R)
  z(2) = -(1/L) * y(1);    // IL' = -(1/L) * UC
  endfunction

  // données :
  t0    =  0;   // temps initial
  tFin  = 0.1;   // temps final
  y0    = [12; 0];  // UC(0) = 12 [V] ; IL(0) = 0 [m/s]], vecteur COLONNE
  N  = 500;      // nombre de points d'évaluation
  at = linspace(t0, tFin, N+1);  // vecteur LIGNE

  // résolution numérique en appelant la fonction  ode  de SciLab
  ay = ode( y0,..     // condition initiale
            t0,..     // x_initial,
            at,..     // série des t(j) pour lesquels  y(j)  sera calculé 
            1e-4,..   // tolérance d'erreur relative
            1e-5,..   // 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 l'erreur absolue et l'erreur relative.

  // graphique de la solution numérique
  win0 = scf(1);
  clf();
  plot(at, ay(1,:), 'b');  // 'b' = blue, pour les tensions aux bornes du condensateur
  plot(at, 1000*ay(2,:), 'r');  // 'r' = red , pour les intensités à travers la self
  axes0 = win0.children;
  title("Tension [V] en bleue et Courant [mA] en rouge", 'fontsize', 4);
  xlabel(axes0, "Temps [s]", 'fontsize', 3, "color", [0, 0, 0]);
  ylabel(axes0, "Tension [V] et Courant [mA]", 'fontsize', 3, "color", [0, 0, 1]);
endfunction
// ex4

function ex5(decal)
//=================
  // Simulation d'un oscillateur  LC, entretenu
  // decal est le décalage de la fréquence d'entretien, relativement
  // a la fréquence de résonnance.
  // decal = 1  => fréquence de résonnance.
  // decal = 1 + 0.1 => fréquence 10% plus haut que le fréquence de résonnance.
  // decal = 1 - 0.15 => fréquence 15% plus bas que la fréquence de résonnance.

  // Définition de l'équation différentielle :  
  // Simulation d'un circuit R L C, c.f. ex. 1 de la série 7
  function z=f(x, y)
  //================
  // x représente le temps
  // y(1)  représente la tension aux bornes du condensateur    UC  en volts
  // y(2)  représente l'intensité du courant traversant la self  IL  en ampères
  // z(1)  représente la dérivée de y(1)
  // z(2)  représente la dérivée de y(2)
  C = 1.0e-6;  // valeur de la capacité en Farads
  L = 4;  // valeur de la self (inductance) en Henri
  R = 15e3;   // valeur de la résistance, en Ohms 
  w0 = 1 / sqrt(L*C);  // pulsation de résonnance  [1/s]
  w = decal*w0;  // décalage par rapport à la pulsation de résonnance.
  // fréquence de résonnance = 1 / (2 * %pi * sqrt(L*C))  [Hz]
  // Période correspondante = 2 * %pi * sqrt(L*C)  [s]
  U0 = 5*cos(w*x);  // générateur de tension
  z(1) = (1/C) * (y(2) - y(1)/R);  
  z(2) = (1/L) * (U0 - y(1));  
  endfunction

  // données :
  t0    =  0;   // temps initial
  tFin  = 0.15;   // temps final
  y0    = [0; 0];  // UC(0) = 12 [V] ; IL(0) = 0 [m/s]], vecteur COLONNE
  N  = 500;      // nombre de points d'évaluation
  ax = linspace(t0, tFin, N+1);  // vecteur LIGNE

  // résolution numérique en appelant la fonction  ode  de SciLab
  ay = ode( y0,..     // condition initiale
            t0,..     // x_initial,
            ax,..     // série des x(j) pour lesquels  y(j)  sera calculé 
            1e-4,..   // tolérance d'erreur relative
            1e-5,..   // 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 l'erreur absolue et l'erreur relative.

  // graphique de la solution numérique
  win0 = scf(2);
  clf();
  plot(ax, ay(1,:), 'b');  // 'b' = blue, pour les tensions aux bornes du condensateur
  plot(ax, 1000*ay(2,:), 'r');  // 'r' = red , pour les intensités à travers la self
  axes0 = win0.children;
  title("Tension [V] en bleue et Courant [mA] en rouge", 'fontsize', 4);
  xlabel(axes0, "Temps [s]", 'fontsize', 3, "color", [0, 0, 0]);
  ylabel(axes0, "Tension [V] et Courant [mA]", 'fontsize', 3, "color", [0, 0, 1]);
endfunction
// ex5
