Discussion :

[rostomus]

Il faut d'abord savoir que les fonctions ode résolvent les EDs de premier ordre, et aussi qu'elles n'expriment pas les résultats sous forme d'expressions mathématique; elles renvoient un vecteur colonne représentant la variable "T" (temps en générale) et une matrice "Y" dont les colonnes sont les solutions.

La fonction ode45 est la première fonction à essayer .

Pour résoudre un système d'équations différentielles :

• vous devez créer une fonction m.file qui définit le système, la fonction a en entrée un scalaire "t" et un vecteur "y" et en sortie un vecteur "dy" qui représente la dérivée.
  Par exemple si on veut résoudre le système:
  
  
  
  
  On crée la fonction suivante:
  

  Code :

  Sélectionner tout - Visualiser dans une fenêtre à part

  1
  2
  3
  4
  5
  function dy = rigid(t,y)
  dy = zeros(3,1);    % a column vector
  dy(1) = y(2) * y(3);
  dy(2) = -y(1) * y(3);
  dy(3) = -0.51 * y(1) * y(2);

• Ensuite, on utilise la fonction odes45 comme suit:
  

  Code :

  Sélectionner tout - Visualiser dans une fenêtre à part

  [T,Y] = ode45(@odefun,plage_t,val_init,options);

  avec "odefun" est le nom de la fonction créée, "plage_t" est un vecteur représente l’intervalle de "t", et "val_init" est un vecteur donnant les valeurs initiales des solutions. pour notre exemple :
  

  Code :

  Sélectionner tout - Visualiser dans une fenêtre à part

  1
  2
  [T,Y] = ode45(@rigid,[0 12],[0 1 1]);
  plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')

  ce qui donne après l'affichage :
  
  
  

Pour plus d'information sur les fonction ode, voir ici

[Invité]

Les fonctions ode ne résolvant que des équations différentielles d'ordre 1, il est tout de même possible de résoudre des EDs d'ordres supérieurs en décomposant le système de sorte à n'avoir que des équations d'ordre 1.

Prenons le système suivant (a étant différent de 0, sinon on se ramène à une équation du premier degré) :

En posant : , nous avons : et

Vient alors le système :

On pourra donc écrire de la même façon la fonction « odefun » comme suit :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
function dy = odefun(t,y)
a = ...
b = ...
c = ...
d = ...
 
dy(1,1) = (d - b*y(1) - c*y(2))/a;
dy(2,1) = y(1);

et l'utiliser :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
plage_t = ...
yp0 = ...
y0 = ...
 
[T,Y] = ode45(@odefun, plage_t, [yp0 y0]);
% Y(:,1) <=> y'(t)
% Y(:,2) <=> y(t)
plot(T, Y(:,1), T, Y(:,2), 'r')
legend('y\prime(t)', 'y(t)')

Cette méthode est bien sûr extensible à des ordres plus élevés.

[maddy22]

Merci beaucoup pour ce tutoriel super clair !
J'aurais cependant une question:
comment fait-on si dans un système on a des équations du type : y1'+ay2'=by2+cy1 par exemple?