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

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

// Exercice 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  
function ex1()
//============
av(1) = 0.0; // Vitesse initiale
at(1) = 0.0; // temps initial
dt = 0.1;  // pas de la méthode numérique de résolution d'équ. diff.
i = 1;
while (at(i) < 10) do // tant que 10 secondes n'ont pas été atteint, calcule ...
   at(i+1) = i * dt;
   av(i+1) = av(i) + dt * (9.81 - (1.5/3)*av(i));  // évolution de la vitesse
   i = i + 1;
end;

scf(1);
clf();
plot(at, av);
endfunction


// Exercice 2.  Équation différentielle, évolution du volume d'eau d'un bassin
// Équation différentielle d'évolution de la hauteur.:
// h'(t) = 2 - h/t) / 6 ,  h(0) = 0  
// 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 ex2()
//============
ah(1) = 0.0; // hauteur initiale
at(1) = 0.0; // temps initial
dt = 0.1;  // pas de la méthode numérique de résolution d'équ. diff.
i = 1;
while (ah(i) < 11.5) do // tant que la hauteur est assez basse
   at(i+1) = i * dt;
   ah(i+1) = ah(i) + dt * ( 2 - ah(i) / 6); // hauteur au temps at(i+1)
   i = i + 1;
end;

scf(2);
clf();
plot(at, ah); // graphique d'évolution de la hauteur en fonction du temps.
disp(at($), ah($));  // Il faut environ 19 secondes pour atteindre une hauteur de 11.5 dm.
endfunction


// Exercice 3.  Évolution d'un oscillateur
// Décrit par une équation différentielle du second ordre :
// x(t)'' = -k * x(t)
// ou par un système de deux équations différentielles du première ordre :
// x'(t) = v(t)
// v'(t) = -k * x(t)
function ex3()
//============
av(1) = 0.0; // Vitesse initiale
ax(1) = 0.1;   // Position initiale
k = 2.0;     // constante de rappel
at(1) = 0.0; // temps initial
T = 20.0; // temps en secondes
N = 1000; // Nombre de pas
dt = T / N;  // pas de la méthode numérique de résolution d'équ. diff.
for i=1:N
   at(i+1) = i * dt;
   ax(i+1) = ax(i) + dt * av(i);  // évolution de la position
   av(i+1) = av(i) - dt * k * ax(i); // évolution de la vitesse
end;

scf(3);
clf();
plot(at, ax);
endfunction

function ex3b()
//============
// Même fonction que la  ex3(),  écrite sous forme vectorielle
axv(:, 1) = [0.1, 0.0]; // Position initiale et Vitesse initiale
k = 2.0;     // constante de rappel
at(1) = 0.0; // temps initial
T = 20.0; // temps en secondes
N = 1000; // Nombre de pas
dt = T / N;  // pas de la méthode numérique de résolution d'équ. diff.
for i=1:N
   at(i+1) = i * dt;
   axv(:, i+1) = axv(:, i) + dt * [axv(2, i) ; -k*axv(1, i)];
end;

scf(3);
plot(at', axv(1, :));   // le  '  transforme le vecteur colonne en ligne
endfunction


// Exercice 4.  Évolution d'un oscillateur avec frottement
// Décrit par une équation différentielle du second ordre :
// x(t)'' = -k * x(t) -c * v(t))
// ou par un système de deux équations différentielles du première ordre :
// x'(t) = v(t)
// v'(t) = -k * x(t) -c * v(t)
function ex4()
//============
axv(:, 1) = [0.1, 0.0]; // Position initiale et Vitesse initiale
k = 2.0;   // constante de rappel
c = 0.3;  // Coef. de frottement proportionnelle à la vitesse
at(1) = 0.0; // temps initial
T = 20.0; // temps en secondes
N = 2000; // Nombre de pas
dt = T / N;  // pas de la méthode numérique de résolution d'équ. diff.
for i=1:N
   at(i+1) = i * dt;
   axv(:, i+1) = axv(:, i) + dt * [axv(2, i) ; -k*axv(1, i)-c*axv(2, i)];
end;

scf(4);
clf();
plot(at', axv(1, :));   // le  '  transforme le vecteur colonne en ligne
endfunction


