// ODE_09_o4_d4_Oscillateurs_couples.sce
// Résolution d'une équation différentielle
// par la méthode de Runge-Kutta qui est d'ordre 4.
// Deux oscillateurs représentés par une masse oscillante
// au bout d'un ressort, sont couplés par un autre ressort
// reliant les deux masses l'une à l'autre.

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 de la masse 1
// y(2) = vitesse  de la masse 1
// y(3) = position de la masse 2
// y(4) = vitesse  de la masse 2
// z(j) = y'(j)
// z(1)= = dérivée de la position de la masse 1
// z(2)= = dérivée de la vitesse  de la masse 1
// z(3)= = dérivée de la position de la masse 2
// z(4)= = dérivée de la vitesse  de la masse 2
m1 = 3;            // masse en [kg] de l'objet 1 
m2 = 3;            // masse en [kg] de l'objet 2 
k01 = 10;  // coefficient de rappel en [N/m] de l'objet 1
k02 = 10;  // coefficient de rappel en [N/m] de l'objet 2
k12 =  1;  // coefficient de rappel du ressort entre deux.
kFrot = 0.00; // coefficient de frottement prop. à la vitesse en [N*s/m]
z(1) = y(2);
z(3) = y(4);
z(2) = -(k01/m1) * (y(1) - 0.10) + (k12/m1) * (y(3) - y(1) - 0.10) - kFrot*y(2)/m1;
z(4) = -(k02/m2) * (y(3) - 0.20) - (k12/m2) * (y(3) - y(1) - 0.10) - kFrot*y(4)/m2;
endfunction

function y=Ordre4(x0,xFin,y0,nMax)
//================================
// pour la condition initiale  y(x0) = y0
// résoud l'équation différentielle y'(x)=f(x,y(x))
// pour  x  variant de  x0  à  xFin
// retourne un vecteur de  nMax+  valeurs de  y.
y    = zeros(length(y0), nMax+1);
y(:,1) = y0;
dx    = (xFin - x0) / nMax;

for j=1:nMax do
  x = x0 + (j-1)*dx
  k1 = f(x      , y(:, j));
  k2 = f(x + dx/2, y(:, j) + k1 * dx/2);
  k3 = f(x + dx/2, y(:, j) + k2 * dx/2);
  k4 = f(x + dx  , y(:, j) + k3 * dx  );
  y(:, j+1) = y(:, j) + dx*(k1 + 2*k2 + 2*k3 + k4)/6;
  end;
endfunction

// Utilisation de la résolution
// données :
a  = 0;
b  = 36.1;
y0 = [0.14; 0; 0.20-0.00; 0]; // pos. 1, vitesse 1, pos. 2, vitesse 2
N  = 2000;
h = (b-a)/N;
x = a:h:b;

// résolution numérique
y  = Ordre4(a,b,y0,N);

// 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
          );
// 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
clf;
plot(x, y(1,:)); // position de l'objet 1 en fonction du temps
plot(x, y(3,:)); // position de l'objet 2 en fonction du temps
disp( string(y(1,N+1)) + "  " + string(y(2,N+1)) + "  " +..
      string(y(3,N+1)) + "  " + string(y(4,N+1)));
      
// poour tester la théorie avec k12 = 0      
//yex = 0.1 + 0.04 * cos(sqrt(10/3)*x);
//plot(x, yex, "g");

//plot(x, y(3,:)-y(1,:), "r");
//plot(x, y(3,:)+y(1,:), "r");
 
return;
plot(x, y(3,:) - y(1,:), "g"); // distance entre les deux objets
plot(x, 0.5*10*(y(1,:)-0.10)^2 + 0.5*3*y(2,:)^2); // énergie de l'objet 1
plot(x, 0.5*10*(y(3,:)-0.20)^2 + 0.5*3*y(4,:)^2); // énergie de l'objet 2
plot(x, 0.5*10*(y(1,:)-0.10)^2 + 0.5*3*y(2,:)^2 +..
        0.5*10*(y(3,:)-0.20)^2 + 0.5*3*y(4,:)^2, "r"); // énergie de l'objet 1 + objet 2
        
