Résolution des systèmes d'équations différentielles ordinaires avec Matlab


Pour résoudre des équations différentielles avec Matlab, il faut procéder en quatre étapes. Celles-ci sont illustrées dans le cas de la simulation de l'amortissement sec d'un oscillateur (§ 3.9.2 du cours).

1. Mise sous forme canonique de l'équation différentielle.

Matlab est capable de résoudre numériquement des systèmes d'équations différentielles du première ordre. Pour résoudre les équations de la mécanique, lesquelles sont généralement d'ordre deux, il est donc nécessaire de transformer chacune des équations d'ordre deux en deux équations du premier ordre. À chaque équation du type

x''(t) = f (t,x(t),x'(t))

on substituera donc le système équivalent

x'(t) = y(t)

y'(t) = f(t, x(t), y(t))

dont la nouvelle inconnue y(t) s'identifie à la vitesse x'(t).

Après transformation des équations différentielles, le système peut généralement s'écrire sous la forme canonique

X'(t) = F (t, X(t) )

où X(t) désigne la matrice colonne des inconnues et F représente la matrice colonne des seconds membres.

 Cas de l'oscillateur amorti par frottement sec :

 m x''(t) + k x(t) = - mu m g sign(x'(t))

(si x'<>0 et | k x | > mu m g)

 devient

 x'(t) = y(t)

y'(t) = - mu g sign(x'(t)) - k x(t)/ m

2. Programmation du second membre.

Pour permettre la résolution de ce système en utilisant les outils de Matlab, il faut définir une fonction implémentant le calcul des seconds
membres du système sous la forme d'une matrice colonne. Cette fonction est généralement décrite dans un fichier séparé. Ses deux premiers arguments doivent impérativement être la variable indépendante (le temps) et la matrice colonne des variables dépendantes (les inconnues X). Des arguments supplémentaires peuvent être utilisés pour passer des paramètres. La valeur de la fonction doit être la matrice colonne des seconds membres de l'équation différentielle (la matrice colonne F de la forme canonique).

 Cas de l'oscillateur amorti par frottement sec :

function xdot=fs(t,x,flag,om,g,mu)
% Fonction correspondant à Coulomb.m

xdot(1) = x(2);
xdot(2) = - om^2*x(1)-mu*g*sign(x(2));

if ((abs(x(2)) <= 0.00000001) & (abs(om^2*x(1))<= mu*g) )
xdot(2)=0;
end

xdot=xdot';

 xdot désigne la valeur de la fonction
La fonction utilise les paramètres om, g et mu

Première équation
Seconde équation

Critère d'arrêt correspondant à la condition
x'<>0 et | k x | > mu m g

Transposition de la sortie pour obtenir une matrice colonne

3. Résolution du système d'équations.

La résolution proprement dite de l'équation différentielle est réalisée au moyen d'une instruction du type

[t,x]=ode23('nom_du_fichier_décrivant_la_fonction',tspan,x0);

où tspan désigne le domaine d'intégration et x0 représente la condition initiale. Le résultat se présente sous la forme de colonnes t et x qui reprennent respectivement les différents instants auxquels la solution a été calculée et les valeurs correspondantes des variables dépendantes.

Différents solveurs (autres que ode23) sont disponibles pour traiter des problèmes difficiles.

Des paramètres peuvent être utilisés en complétant l'instruction d'appel au solveur :

[t,x]=ode23('nom_du_fichier_décrivant_la_fonction',tspan,x0,[], paramètre 1, paramètre 2, ...);

 Cas de l'oscillateur amorti par frottement sec :

% Liste des paramètres
k = 0.1 ;
g = 9.81 ;
m = 0.05 ;
mu = 0.01 ;
om = sqrt(k/m) ;
t0=0; tf=3*2*pi/om;
tspan = [t0 tf];
x0 = [1 0]';

% Résolution de l'équation différentielle - frottement sec

[t,x]=ode23('Coulomb_func',tspan,x0,[],om,g,mu);

Initialisation des valeurs des paramètres

 

 

Définition de la durée de l'intégration
Conditions initiales

Coulomb_func est ici le nom du fichier décrivant le second membre de l'équation (Tableau précédent)


4. Représentation graphique des résultats

Les résultats sont obtenus sous la forme des matrices colonnes t et x. Leur représentation graphique est possible au moyen de la commande plot.

 Cas de l'oscillateur amorti par frottement sec :

xlim([t0 tf]);

plot (t,x(:,1))

xlabel('t'); ylabel('x');


legend(strcat('\mu=' ,num2str(mu)));

Définition des limites de la fenêtre d'affichage (t0 <= t <= tf)

Représentation graphique

Définition des labels portés par les axes

Affichage d'une légende

 

P.S. Les fichiers matlab correspondant à cet exemple (Coulomb.m et Coulomb_func.m) sont disponibles sur MathServer dans le répertoire /u/MC/matlab. Pour l'execution, il suffit de taper la commande Coulomb à l'invite de commande de Matlab. Remarquez que ces fichiers ne peuvent être modifiés tels quels. Ils doivent être copiés au préalable dans votre répertoire.


Travaux personnels :

L'approche décrite ci-dessus peut être appliquée pour résoudre numériquement les équations différentielles rencontrées dans les applications.