// Série 07,  résolution numérique d'équations différentielles
// de dimensions supérieur à 1

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 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
   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(1);
   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
// ex1

function ex2(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
// ex2


function ex3()
//============
// Résolution d'une équation différentielle
// par la méthode du système.
// Problème de la trajectoire d'une comète tombant sur la Terre.

function z=f(x,y)
//===============
// pour résoudre vectoriellement :
// v' = -(G*MTerre/(d^3)) * [x, z] - (k_frot/m) * ||v|| * v.
// d = sqrt(x^2 + z^2) = la distance au centre de la Terre.
// masse = m = 350 kg
// y(1) = position  x  de l'objet
// y(2) = vitesse   x  de l'objet
// y(3) = position  z  de l'objet
// y(4) = vitesse   z  de l'objet
// La position  y  de l'objet étant fixe, on n'en tient pas compte.
// z(j) = y'(j)
// z(1)= = dérivée de la position en x
// z(2)= = dérivée de la vitesse en x
// z(3)= = dérivée de la position en z
// z(4)= = dérivée de la vitesse en z
G = 6.6726e-11;     // en N*m^2/kg^2
m = 350;            // masse de l'objet en [kg] 
MTerre = 5.974e24;  // masse de la Terre en kg
z(1) = y(2);
z(3) = y(4);
scal = -G*MTerre / sqrt( y(1)^2 + y(3)^2 )^3;
vitesse   = sqrt( y(2)^2 + y(4)^2 );  // vitesse de l'objet
z(2) = scal * y(1) -(3e-7/m) * vitesse * y(2);
z(4) = scal * y(3) -(3e-7/m) * vitesse * y(4);
endfunction

// Utilisation de la résolution
// données :
RTerre = 6.371e6; // Rayon de la Terre en [m]
a  = 0;
b  = 33000;
y0 = [0; 3000; 3*RTerre; 0]; // pos. x, vitesse x, pos. z, vitesse z
N  = 1000;
h = (b-a)/N;
x = a:h:b;

// résolution numérique
y  = ode( y0,..    // condition initiale
          a,..     // x_initial,
          x,..     // 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
          );
// Le maximum entre l'erreur absolue et l'erreure relative sera recherché,
// la le minimum entre ces deux données.

// graphique de la solution numérique
win0 = scf(3);
clf();
isoview(-1.5, 3, -1.5, 3);

plot(y(1,:)/RTerre, y(3,:)/RTerre);  // trajectoire de la comète
axes0 = win0.children;
title("Trajectoire de la comète", 'fontsize', 4);
xlabel(axes0, "Position X [RTerre]", 'fontsize', 3, "color", [0, 0, 1]);
ylabel(axes0, "Position Y [RTerre]", 'fontsize', 3, "color", [0, 0, 1]);
endfunction
// ex3


function ex4()
//============
// Résolution d'une équation différentielle
// par la méthode du système.
// Problème de deux corps qui se trouvent dans l'espace,
// chacun attire l'autre selon la loi de la gravitation générale.
// Ils sont isolé de toute autre interaction.

function z=f(x,y)
//===============
// pour résoudre vectoriellement :
// v1' = -(k01/m1) * (x1 - 0.10) + (k12/m1) * (x2 - x1 - 0.01)
// v2' = -(k02/m2) * (x2 - 0.20) - (k12/m2) * (x2 - x1 - 0.01)
// La position au repos de l'objet 1 est à 0.10 [m]
// La position au repos de l'objet 2 est à 0.20 [m]
// y(1) = position x de la masse 1
// y(2) = vitesse  x de la masse 1
// y(3) = position y de la masse 1
// y(4) = vitesse  y de la masse 1
// y(5) = position x de la masse 2
// y(6) = vitesse  x de la masse 2
// y(7) = position x de la masse 2
// y(8) = vitesse  x de la masse 2
// z(j) = y'(j)
G = 6.6726e-11;   // en N*m^2/kg^2
m1 = 3e14;        // masse en [kg] de l'objet 1 
m2 = 3e14;        // masse en [kg] de l'objet 2 
Force = G * m1 * m2 / sqrt( (y(5) - y(1))^2 + (y(7) - y(3))^2 )^3;
z(1) = y(2);
z(3) = y(4);
z(5) = y(6);
z(7) = y(8);
z(2) = Force/m1 * ( y(5) - y(1) );
z(4) = Force/m1 * ( y(7) - y(3) );
z(6) = Force/m2 * ( y(1) - y(5) );
z(8) = Force/m2 * ( y(3) - y(7) );
endfunction

// Utilisation de la résolution
// données :
a  = 0;
b  = 125;
y0 = [-100; 1; 0; 6;..  // pos_x. 1, vitesse_x 1, pos_y. 1, vitesse_y 1
       100; 0; 0; -6]; // pos_x. 2, vitesse_x 2, pos_y. 2, vitesse_y 2
N  = 400;
N  = 1000;
h = (b-a)/N;
x = a:h:b;

// résolution numérique
y  = ode( y0,..    // condition initiale
          a,..     // x_initial,
          x,..     // 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
          );
// Le maximum entre l'erreur absolue et l'erreure relative sera recherché,
// la le minimum entre ces deux données.

// graphique de la solution numérique

win0 = scf(4);
clf();
plot(y(1,:), y(3,:));  // trajectoire du corps 1
plot(y(5,:), y(7,:), "r");  // trajectoire du corps 2
axes0 = win0.children;
//axes0.isoview = "on";
endfunction
// ex4

