// Série 04,  Simulations numériques de problèmes physique

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
ax(1) = 0.0; // Position initiale
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
   ax(i+1) = ax(i) + dt * av(i);  // évolution de la position
   i = i + 1;
end;

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

scf(2);
clf();
plot(at, ax);
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 amorti
// ou par un système de deux équations différentielles du première ordre :
// x ' (t) = v(t)
// v ' (t) = -k/m * x(t) - c/m * v(t)
function ex3()
//============
at(1) = 0.0; // temps initial
ax(1) = 0.1;   // Position initiale
av(1) = 0.0; // Vitesse initiale
m = 0.80;  // masse du corps
k = 2.0;     // constante de rappel
c = 0.3;  // Coef. de frottement proportionnelle à la vitesse
Tmax = 20.0; // temps en secondes
dt = 0.02;  // pas de la méthode numérique de résolution d'équ. diff.
t = 0;  // temps
i = 1;  // indice pour mémoriser l'évolution
while (t < Tmax)
   t = i * dt;
   at(i+1) = t;
   a_de_t = -(k * ax(i) + c * av(i)) / m;
   ax(i+1) = ax(i) + dt * av(i);  // évolution de la position
   av(i+1) = av(i) + dt * a_de_t; // évolution de la vitesse
   i = i + 1;
end;

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

//ex3();

// 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/m * x(t) -c/m * v(t)
function ex3b()
//============
// Même fonction que la  ex3(),  écrite sous forme vectorielle
at(1) = 0.0; // temps initial
axv(:, 1) = [0.1, 0.0]; // Position initiale et Vitesse initiale
m = 0.80;  // masse du corps
k = 2.0;   // constante de rappel
c = 0.3;  // Coef. de frottement proportionnelle à la vitesse
Tmax = 20.0; // temps en secondes
dt = 0.02;  // pas de la méthode numérique de résolution d'équ. diff.
t = 0;  // temps
i = 0;  // indice pour mémoriser l'évolution
while (t < Tmax)
   i = i + 1;
   t = i * dt;
   at(i+1) = t;
   a_de_t = -(k * axv(1, i) + c * axv(2, i)) / m;  // accélération au temps (t-dt)
   axv(:, i+1) = axv(:, i) + dt * [ axv(2, i) ;.. // axv[2,i] = vitesse au temps (t-dt)
                                    a_de_t];      // Vecteur colonne, de 2 lignes.
end;

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

//ex3b();

return;
N=100
ss=0
for ii=1:N do
    ss=ss+9^(N-ii)*10^ii
end;
disp( 1- ss/10^(N+1));

ss=0
for ii=0:(N-1) do
    ss=ss+0.9^(ii)
end;
disp( 1 - 0.1*ss);

ss = (1 - 0.9^N) / (1 - 0.9);
disp( 1 - 0.1*ss);
