Équations différentielles

Problème : on cherche à résoudre un système d'équations différentielles ordinaires (EDO, en anglais ODE) du premier ordre :

\[\begin{align*} \dot{y}_1 &= f_1(y_1,y_2,\dots, y_n,t) \\ \dot{y}_2 &= f_2(y_1,y_2,\dots, y_n,t) \\ &\dots \nonumber \\ \dot{y}_n &= f_n(y_1,y_2,\dots, y_n,t) \end{align*} \]

avec les conditions initiales \(y_1(0) = y_{10} \quad \dots \quad y_n(0) = y_{n0}\)

Rappelons au lecteur non averti que tout système ou équation différentielle d'ordre supérieur peut se ramener simplement à cette forme canonique, utilisée dans tous les solveurs d'EDO.

On voit donc que la définition d'un tel système repose sur la définition de \(n\) fonctions de \(n+1\) variables. Ces fonctions devront être programmées dans une fonction MATLAB sous la forme canonique suivante :

function ypoint = f (t, y)
 
   ypoint(1) = une expression de y(1), y(2) ... y(n) et t
   ...
   ypoint(n) = une expression de y(1), y(2) ... y(n) et t
   
   ypoint = ypoint(:);
 
end

On remarquera que les \(y_i\) et les \(\dot y _i\) sont regroupés dans des vecteurs, ce qui fait que la forme de cette fonction est exploitable quel que soit le nombre d'équations du système différentiel.

La dernière ligne est nécessaire ici, car la fonction doit renvoyer un vecteur colonne et non un vecteur ligne.

Évidemment, sachant que les expressions des dérivées doivent être stockées dans un vecteur colonne, on peut écrire directement :

function ypoint = f (t, y)
 
   ypoint(1, 1) = une expression de y(1), y(2) ... y(n) et t
   ...
   ypoint(n, 1) = une expression de y(1), y(2) ... y(n) et t
   
 
end

Ensuite, pour résoudre cette équation différentielle, il faut appeler un solveur et lui transmettre au minimum :

  • le nom de la fonction.

  • les bornes d'intégration ( \(t_{min}\) et \(t_{max}\) ).

  • les conditions initiales.

Le solveur fournit en sortie un vecteur colonne représentant les instants d'intégration \(t\), et une matrice dont la première colonne représente les \(y_1\) calculés à ces instants, la deuxième les \(y_2\), et la \(n^{i\grave{e}me}\) les \(y_n\).

L'appel du solveur prend donc en général la forme suivante :

[t, y] = ode45 (@f, [tmin tmax], [y10 ; y20 ; ... ; yn0] );
y1 = y(:,1);
y2 = y(:,2);
...
yn = y(:,n);
plot(t, y1, t, y2)     % par exemple on trace y1(t) et y2(t)
plot(y1,y2)           % ou bien y2(y1) (plan de phase pour les oscillateurs)

Les lignes y1 = ... servent à extraire les différentes fonctions \(y_i\) dans des colonnes simples.

Nous avons utilisé ici ode45 qui est un Runge-Kutta-Merson imbriqué d'ordre 4 et 5. C'est le plus courant et celui par lequel il faut commencer, mais il en existe d'autres, en particulier ode15s adapté aux systèmes raides (voir la doc).

Les spécialistes s'étonneront de ne pas avoir à spécifier d'erreur maximale admissible, relative ou absolue. Sachez que MATLAB prend une erreur relative max de \(10^{-4}\) par défaut, et qu'il est toujours possible de modifier cette valeur, ainsi que bien d'autres paramètres grâce à la routine de gestion des options odeset.

Exemple

Il est temps de passer à un exemple.

On considère l'équation de Matthieu amortie :

\[\ddot{y} + b\dot{y} + a \left( 1+\epsilon \cos \left( t\right) \right) y = 0\]

où \(a\), \(b\) et \(\epsilon\) sont des paramètres. On prend comme conditions initiales \(y(0) = 10^{-3}\) et \(\dot{y}(0) = 0\). En posant \(y_1 = y\) et \(y_2 = \dot{y}\)

on se ramène à la forme canonique :

\[\begin{align*} \dot{y}_1 &= y_2 \\ \dot{y}_2 &= - b y_2 -a \left( 1+\epsilon \cos \left( t \right) \right) y_1 \end{align*}\]

Écrivons la fonction matthieu définissant cette équation dans un fichier matthieu.m. Dans cet exemple, les paramètres de l'équation devront être passés comme entrées de la fonction :

 function ypoint = matthieu (t, y, a, b, epsilon) 

  ypoint(1,1) = y(2);
  ypoint(2,1) = -b*y(2) -a*(1+epsilon*cos(t))*y(1);

end

Pensez à mettre des ; à la fin de chaque ligne si vous ne voulez pas voir défiler des résultats sans intérêt.

La séquence d'instructions (à mettre dans un autre fichier .m) qui appelle le solveur sera par exemple :

% Paramètres
a = 1;
b = 0.1;
epsilon = 1;

% Temps final
tfinal = 10*pi;

% Conditions initiales
y01 = 1e-3;    
y02 = 0;

[t,y] = ode45( @(t, y) matthieu(t, y, a, b, epsilon), [0 tfinal], [y01 y02]); % Résolution

y1 = y(:,1);  	% Extraction de y1 et y2
y2 = y(:,2);

subplot(221)
plot(t,y1)    	% y1 fonction du temps 
              		%      (représente y(t) pour l'EDO originale)
subplot(222)
plot(t,y2)    	% y2 fonction du temps
              		%      (représente dy(t)/dt pour l'EDO originale)
subplot(212)
plot(y1,y2)   	% Plan de phase

Ici, la syntaxe de fonction anonyme ( @(t, y) matthieu(t, y, a, b, epsilon)) permet de se ramener à une fonction dépendant que de \(t\) et de \(y\), pour les paramètres \(a\), \(b\) et \(epsilon\) fixés.

Ce programme trace la figure suivante qui représente les grandeurs \(y(t)\) et \(\dot y(t)\) de l'équation originale en fonction du temps, plus le plan de phase. Au passage, on retrouve bien l'instabilité des solutions de l'équation de Matthieu pour les valeurs des paramètres choisis.

Résultat obtenu pour l'équation de Matthieu avec ode45
Résultat obtenu pour l'équation de Matthieu avec ode45

Remarque

Il est naturellement possible de définir le système d'équations différentielles à résoudre par l'intermédiaire d'une fonction anonyme et non pas avec une fonction externe.

Avec une fonction anonyme, l'exemple précédent est résolu ainsi :

a=1;
b=0.1;
epsilon=1;
%
fMatthieu= @(t,y) [y(2); -b*y(2)-a*(1+epsilon*cos(t))*y(1)];
[t,y] = ode45(fMatthieu, [0 10*pi], [1e-3 0]);