// Série 05,  résolution numérique d'équations différentielles.
// Écriture de fonctions générales de résolution numérique d'équations différentielles.

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

function y=Euler_o1_d1(x0, xFin, y0, nMax, f)
//===========================================
// o1, signifie que l'ordre numérique est de 1, soit peu précis.
// d1, signifie que la dimension de  y0  et de  y  est de  1, soit assez limité.
// Résolution numérique de l'équation différentielle y'(x) = f(x, y(x))
// y(x0) = y0  est la condition initiale
// retourne un vecteur de  nMax+1  valeurs de  y.
// f  est la fonction décrivant l'équ. diff.
y(1) = y0;
dx = (xFin - x0) / nMax;
for j=1:nMax do
   x = x0 + (j-1) * dx;
   y(j+1) = y(j) + dx * f(x, y(j));
   end;
endfunction


// 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()
//============
   function z=f(x, y)
   //================
   // z  représente  y'
   z = 9.81 - (1.5/3)*y;  // car  y' = f(x, y) = 9.81 - (1.5/3) * y
   endfunction
   
   // Résolution numérique
   t0 = 0;
   tFin = 10;
   N = 100;
   ay = Euler_o1_d1( t0..    // temps initial
                    ,tFin..   // temps final
                    ,0..    // vitesse initiale
                    ,N..  // nombre de points = nMax
                    ,f);    // f  est la fonction définissant l'équ. diff.
   
   // graphique de la solution numérique
   ax = linspace(t0, tFin, N+1)';   // le  '  transforme le vecteur ligne en colonne
   scf(1);
   clf();
   plot(ax, ay, 'k');
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()
//============
   function z=f(x, y)
   //================
   // z  représente  y'
   z = 2 - y/6;  // car  y' = f(x, y) = 2 - y/6
   endfunction
   
   tFin = 20;
   N = 200;
   // Résolution numérique
   ay = Euler_o1_d1( 0..    // temps initial
                    ,tFin..   // temps final
                    ,0..    // hauteur initiale
                    ,N..  // nombre de points = nMax
                    ,f);    // f  est la fonction définissant l'équ. diff.
   
   // graphique de la solution numérique
   ax = linspace(0, tFin, N+1)';   // le  '  transforme le vecteur ligne en colonne
   scf(2);
   clf();
   plot(ax, ay, 'k');
endfunction

//###################################################3

function y=Euler_o1_dn(x0, xFin, y0, nMax, f)
//===========================================
// o1, signifie que l'ordre numérique est de 1, soit peu précis.
// dn, signifie que la dimension de  y0  et de  y  peut être quelconque.
// Résolution numérique de l'équation différentielle y'(x) = f(x, y(x))
// y(x0) = y0  est la condition initiale
// retourne un vecteur de  nMax+1  valeurs de  y.
// f  est la fonction décrivant l'équ. diff.
// y0  et  f  doivent être des vecteurs colonne.
y(:, 1) = y0;
dx = (xFin - x0) / nMax;
for j=1:nMax do
   x = x0 + (j-1) * dx;
   y(:, j+1) = y(:, j) + dx * f(x, y(:, j));
   end;
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()
//============
   function z=f(x, y)
   //================
   // x représente le temps
   // y(1)  représente la position
   // y(2)  représente la vitesse
   // z(1)  représente la dérivée de la position par rapport au temps
   // z(2)  représente la dérivée de la vitesse par rapport au temps
   // y et z sont des vecteurs colonne
   k = 2.0;
   z(1) = y(2);     // car  position ' = vitesse
   z(2) = -k*y(1);  // car  vitesse ' = -k * position
   endfunction
   
   // données :
   t0    = 0;    // temps initial
   tFin  = 20;   // temps final
   y0 = [0.1; 0];  // position ; vitesse initiale, vecteur colonne
   N  = 1000;      // nombre de points d'évaluation
   // Résolution numérique
   ay = Euler_o1_dn( t0, tFin, y0, N, f);
   
   // graphique de la solution numérique
   at = linspace(t0, tFin, N+1);  // est un vecteur ligne
   scf(3);
   clf();
   plot(at, ay(1, :), 'k');
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()
//============
   function z=f(x, y)
   //================
   // x représente le temps
   // y(1)  représente la position
   // y(2)  représente la vitesse
   // z(1)  représente la dérivée de la position par rapport au temps
   // z(2)  représente la dérivée de la vitesse par rapport au temps
   // y et z sont des vecteurs colonne
   k = 2.0;   // constante de rappel
   c = 0.3;  // Coef. de frottement proportionnelle à la vitesse
   z(1) = y(2);     // car  position ' = vitesse
   z(2) = -k*y(1) -c*y(2);  // car  vitesse ' = -k * position - c * vitesse
   endfunction
   
   // données :
   t0    = 0;    // temps initial
   tFin  = 20;   // temps final
   y0 = [0.1; 0];  // position ; vitesse initiale, vecteur colonne
   N  = 1000;      // nombre de points d'évaluation
   // Résolution numérique
   ay = Euler_o1_dn( t0, tFin, y0, N, f);
   
   // graphique de la solution numérique
   at = linspace(t0, tFin, N+1);  // est un vecteur ligne
   scf(4);
   clf();
   plot(at, ay(1, :), 'k');
endfunction


// Exercice 5.  Test de la précision de la méthode d'Euler
// x'(t) =  v(t)
// v'(t) = -x(t)
// x(0) = 0;  v(0) = 1;
// La solution exacte est x(t) = sin(t) ; v(t) = cos(t)
// On peut donc comparer la méthode numérique avec la solution exacte.
function ex5()
//============
   function z=f(x, y)
   //================
   // x représente le temps
   // y(1)  représente la position
   // y(2)  représente la vitesse
   // z(1)  représente la dérivée de la position par rapport au temps
   // z(2)  représente la dérivée de la vitesse par rapport au temps
   // y et z sont des vecteurs colonne
   z(1) = y(2);   // car  position ' =  vitesse
   z(2) = -y(1);  // car  vitesse '  = -position
   endfunction
   
   // données :
   t0    = 0;        // temps initial
   tFin  = 2*%pi;   // temps final
   y0 = [0; 1];  // position ; vitesse initiale, vecteur colonne
   N  = 200;      // nombre de points d'évaluation
   // Résolution numérique
   ay = Euler_o1_dn( t0, tFin, y0, N, f);
   
   // graphique de la solution numérique
   at = linspace(t0, tFin, N+1);  // est un vecteur ligne
   scf(3);
   clf();
   plot(at, ay(2, :), 'k');
   plot(at, cos(at), 'r');
   
   // Affiche la précision de la méthode d'Euler
   disp("N=" + string(N) + "  N*erreur = " + string(N*(ay(2,$) - cos(at($)))));
endfunction

//ex5;


