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

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

function ex1()
//============
// 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(t,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
  k_frot = 2e-7; // [N * s^2 / m^2] Coef. de frottement.
  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) -(k_frot/m) * vitesse * y(2);
  z(4) = scal * y(3) -(k_frot/m) * vitesse * y(4);
  endfunction

// Utilisation de la résolution
// données :
RTerre = 6.371e6; // Rayon de la Terre en [m]
t_fin  = 40000;
y0 = [0; 3000; 3*RTerre; 0]; // pos. x, vitesse x, pos. z, vitesse z
N  = 1000;
dt = t_fin/N;
at = 0:dt:t_fin;

// résolution numérique
y  = ode( y0,..    // condition initiale
          0,..     // x_initial,
          at,..     // série des at(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(1);
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
// ex1()

function ex2()
//============
// 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(t,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 :
t_fin  = 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  = 1000;
dt = t_fin/N;
at = 0:dt:t_fin;

// résolution numérique
y  = ode( y0,..    // condition initiale
          0,..     // x_initial,
          at,..     // 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(2);
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
// ex2())

function ex3()
//============
// Problème de Lotka Volterra, des Lapins et renards
// Problème du système proie-prédateur de Lotka-Volterra.
// c.f. série 7

  function z=f(t,y)
  //===============
  // pour résoudre  y' = (a y(1) - b y(1) y(2) ; -c y(2) + d y(1) y(2))
  // y(1) = nb de lapins   y(2) = nombre de renards
  // z(1) = y'(1)  z(2) = y'(2)
  a = 1;  
  b = 1; 
  c = 1;
  d = 1;
  z(1) = a*y(1) - b*y(1)*y(2);
  z(2) = -c*y(2) + d*y(1)*y(2);
  endfunction

// Utilisation de la résolution
// données :
t_fin  = 50;
y0 = [1; 0.5];
N  = 1000;
dt = t_fin/N;
at = 0:dt:t_fin;

// résolution numérique
y  = ode( y0,..    // condition initiale
          0,..     // x_initial,
          at,..     // 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
scf(3);
clf();
//plot(at, y(2,:));
plot(y(1,:),y(2,:));
endfunction
// ex3())


function ex4()
//============
// 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és de toute autre interaction.

  function z=f(t,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
  // y( 9) = position x de la masse 3
  // y(10) = vitesse  x de la masse 3
  // y(11) = position y de la masse 3
  // y(12) = vitesse  y de la masse 3
  // 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 
  m3 = 3e12;        // masse en [kg] de l'objet 3  Plus léger ! 
  Force12 = G * m1 * m2 / sqrt( (y(5) - y(1))^2 + (y( 7) - y(3))^2 )^3;
  Force13 = G * m1 * m3 / sqrt( (y(9) - y(1))^2 + (y(11) - y(3))^2 )^3;
  Force23 = G * m2 * m3 / sqrt( (y(9) - y(5))^2 + (y(11) - y(7))^2 )^3;
  z( 1) = y( 2);
  z( 3) = y( 4);
  z( 5) = y( 6);
  z( 7) = y( 8);
  z( 9) = y(10);
  z(11) = y(12);
  z( 2) = Force12/m1 * ( y(5) - y( 1) ) + Force13/m1 * ( y( 9) - y( 1) );
  z( 4) = Force12/m1 * ( y(7) - y( 3) ) + Force13/m1 * ( y(11) - y( 3) );
  z( 6) = Force12/m2 * ( y(1) - y( 5) ) + Force23/m2 * ( y( 9) - y( 5) );
  z( 8) = Force12/m2 * ( y(3) - y( 7) ) + Force23/m2 * ( y(11) - y( 7) );
  z(10) = Force13/m3 * ( y(1) - y( 9) ) + Force23/m3 * ( y( 5) - y( 9) );
  z(12) = Force13/m3 * ( y(3) - y(11) ) + Force23/m3 * ( y( 7) - y(11) );
  endfunction

// Utilisation de la résolution
// données :

// évolution très intéressante !
t_fin  = 2000;  
y0 = [  -50; 0; 0;   9.00;.. // pos_x. 1, vitesse_x 1, pos_y. 1, vitesse_y 1
         50; 0; 0;  -9.15;.. // pos_x. 2, vitesse_x 2, pos_y. 2, vitesse_y 2
       -250; 0; 0;  12];  // pos_x. 3, vitesse_x 3, pos_y. 3, vitesse_y 3
       
t_fin  = 2000;  
y0 = [  -50; 0; 0;   9.00;.. // pos_x. 1, vitesse_x 1, pos_y. 1, vitesse_y 1
         50; 0; 0;  -9.15;.. // pos_x. 2, vitesse_x 2, pos_y. 2, vitesse_y 2
       -250; 0; 0;  12];  // pos_x. 3, vitesse_x 3, pos_y. 3, vitesse_y 3
       
N  = t_fin/2;
dt = t_fin/N;
at = 0:dt:t_fin;

// résolution numérique
y  = ode( y0,..    // condition initiale
          0,..     // x_initial,
          at,..     // série des x(j) pour lesquels  y(j)  sera calculé 
          1e-5,.. // tolérance d'erreur relative
          1e-5,.. // tolérance d'erreur absolue
          f..      // fonction définissant l'équation différentielle
          );

// Animation des deux corps
//============================
scf(3);
clf;
aaa = gca();
aaa.data_bounds=[-700, -700 ; 700, 700];
aaa.isoview="on";  // axes x et y ont la même échelle

plot(y(1,1:2), y(3,1:2));
ee1 = gce();
ee1.children(1).thickness = 9;
plot(y(5,1:2), y(7,1:2), "r");
ee2 = gce();
ee2.children(1).thickness = 9;
plot(y(9,1:2), y(11,1:2), "k");
ee3 = gce();
ee3.children(1).thickness = 9;

//aaa.axes_visible="off";
aaa.margins = [0,0,0,0];

realtimeinit(20/N); // définit un intervalle de temps => temps total = ? [s]
for j=1:N
  realtime(j); // attend le temps définit dans realtimeinit(...)
  drawlater();
  ee1.children(1).data = [y( 1,j), y( 3,j) ; y(1,j)+15, y( 3,j)];
  ee2.children(1).data = [y( 5,j), y( 7,j) ; y(5,j)+15, y( 7,j)];
  ee3.children(1).data = [y( 9,j), y(11,j) ; y(9,j)+15, y(11,j)];
  drawnow();
  
end;

// pour afficher les trois trajectoires
plot(y(1,:), y( 3,:), "b");
plot(y(5,:), y( 7,:), "r");
plot(y(9,:), y(11,:), "k");
endfunction
// ex4()

