Tuto Matlab - ODE 45

//!\\ La compréhension de ce mini-tuto nécessite une connaissance minimale de Matlab //!\\

Salut à toi très cher  amateur (ou amatrice) de simulation numérique, dans ce mini-tuto vous découvrirez comment manipuler une des fonctions les plus utilisées de Matlab, j'ai cité ode45.

Présentation de la fonction ode45

En consultant le help de matlab on découvre que la fonction ode45 reçoit en entrée plusieurs paramètres (4) , et en retourne 2.

 

Les entrées :

  1. Intervalle d'intégration : dans un jargon plus cru ce sont les bornes d'intégrations, A --> B, dans une équation en fonction du temps on va dire que c'est T(début) et T(final).
  2. Équation différentielle : C'est l'équation que ode45 va devoir résoudre. Vous ne devez pas vous tracasser quant à la méthode de résolution à utiliser, Matlab s'en charge pour vous (merci Matlaaaab!!!!!)
  3. Conditions Initiales : vous verrez dans votre cours d'analyse, qu'il sera nécessaire d'avoir des informations sur les conditions de départ pour la résolution de l'équation. Ces conditions dépendent de l'énoncé de votre problème.

    Ex.: si c'est en électricité ce seront les V(0) et I(0), si c'est en mécanique ce sera plutôt les positions et vitesses d'origine, ...
  4. Option : il permet de spécifier quelles sont les conditions à respecter pour effectuer l'intégration. Généralement ce les conditions d'arrêts d'intégration (utile pour détecter un phénomène), la précision, etc ...

Les sorties :

  1. La fonction.
  2. La valeur des couples x,y associés.

Exemple par la pratique:

Pour illustrer ce schéma prenons l'exemple d'une équation de mouvement



Une balle lâché dans l'air à une vitesse nul , hauteur H de 15m

Combien de temps lui faudra t'elle pour atteindre le sol (ou s'écraser)?



Généralement dans l'équation de mouvements on vérifie les forces (ou sollicitation pour un circuit électrique) en jeu. Ici ça sera seulement la gravité (Non monsieur, il n'y a pas de frottement d'air ici, nous sommes dans le vide et puis voilà!).

On suppose que la vitesse initiale est nulle (la balle est lâchée => en t0, v(t0) = 0). La seule force agissant sur le système (mouvement) est la gravité  (vous prendrez l'habitude de dire système, cela fait plus sérieux dans une conversation).

la loi de mouvement régissant le système est très simple :

`frac{d^{2}x}{dt^{2}}` , il s'agit du terme de l'accélération, engendré par la force de gravité (Haelterman dira "le corps soumis au champ gravitationnel G").

La seule équation que ode45 doit intégrer sera donc celle ci:

`frac{d^{2}x}{dt^{2}} = -g`

Nous allons donc coder tout cela dans Matlab (ou octave ou python) afin de résoudre notre problème.

Il faut créer 3 fichiers (il y a moyen de ne le faire qu'en 2 fichiers mais c'est moins joli et surtout moins lisible).

  • Le premier appellera ode45 et servira à l'intégration proprement dite.
  • Le second contiendra l'équation de mouvement 
  • le troisième servira à préciser les conditions d'arrêt



Observons donc ces quelques lignes de code:

calcul.m

x0 = 15

v0 = 0

ConditionInitale = [x0 v0]

T0 = 0

T1 = 2 % en secondes



options = odeset('NonNegative',1); 

[Temps,Solution] = ode45(@equation_mouvement, T0 :0.01 : T1, ConditionInitale,options); % ici ode45 renvoie le resultat dans une structure ce qui permet un traitement facile des résultats

@equation_mouvement indique que nous allons récupérer l'équation du problème dans un autre fichier, nommé... equation_mouvement.m.

T0 et T1 représentent donc l'intervalle d'intégration. Prenez bien garde aux unités de vos variables car Matlab ne détectera pas les éventuelles incohérences.

Notre variable options contient les conditions de résolution du système. Dans ce cas précis nous interdisons au premier paramètre d'être négatif car notre balle ne peut pas s'enfoncer dans le sol.

Nous constatons que ode45 nous ressort un vecteur. Il est donc primordial de bien mettre un vecteur pour récupérer la solution (désolé d'insister aussi lourdement mais c'est malheureusement une faute assez fréquente au début).

Le fichier equation_mouvement.m

function dx = equation_mouvement(t,x)

g = 9.81 %constante de gravitation



% définisson un vecteur de travail pour dire à ode45 où travailler

dx = zero(2,1);



% on initialise dans notre vecteur de travail les conditions initials

dx(1) = x(2)



% L'equation du mouvement tout simplement

dx(2) = -g

 

t représente l'intervalle de temps alors que x est un vecteur contenant les conditions initiales (vitesse et position).

Remarques :
  • nous n'avons pas redéfini x(1) car la vitesse initiale est nulle, et comme notre vecteur a été initialisé à zéro tout est parfait.
  • Contrairement à beaucoup de langages de programmation, Matlab commence à compter à partir de 1 et non de 0. La première case d'un vecteur est donc vecteur(1) et non vecteur(0) (là aussi j'insiste lourdement mais c'est aussi une faute fréquente, surtout auprès de ceux qui ont déjà programmé).

Le fichier arret.m

% cf help de Matlab pour la création de ce fichier. Dans ce cas)ci nous n'avons pas de condition % d'arrêt.



Et voilà, il ne reste maintenant plus qu'à faire exécuter ces quelques lignes de code et la réponse sortira dans la matrice [Temps,Solution] que vous pourrez afficher, mettre en graphique (commande plot), ...

 

Tuto by Momo & Oli

Tags: 

Dernières Publications

Articles Populaires

Nos Coordonnées

Bureau des Étudiants de Polytechnique ASBL (BEP ASBL)

UB1.149 (Campus Solbosch, bâtiment U, porte B, étage 1)

Avenue F.D. Roosevelt, 50 CP 165/02 1050 Bruxelles

bep@bepolytech.be

Tél./Fax : (00-32-2) 650.29.06