SciPy : résoudre une équation qu'on ne sait pas résoudre à la main

Les figures s'affichent dans la page. NumPy, Matplotlib et SciPy sont installés dans l'éditeur du navigateur : l'image apparaît sous le code après exécution. On peut aussi recopier ces exemples sur la machine, ils y fonctionnent à l'identique.

SciPy et le module integrate

C'est quoi SciPy ?

Alors que NumPy fournit les "briques" (les tableaux), SciPy (Scientific Python) fournit les "outils" pour l'ingénieur : traitement du signal, optimisation, statistiques et, ce qui importe ici, l'intégration d'équations différentielles.

L'outil utilisé ici est odeint (Ordinary Differential Equation Integrator). C'est le standard industriel pour résoudre des systèmes dynamiques.

Comment raisonne odeint

Pour utiliser odeint, on doit comprendre que l'ordinateur ne connaît pas la physique. Il ne connaît pas "la gravité" ou "les frottements". Il ne connaît qu'une seule chose : La Pente (la dérivée).

La résolution numérique fonctionne toujours selon le même schéma en 3 étapes :

• L'État actuel (S) : "Je suis ici, à cette vitesse."
• La Dérivée (S′) : "Selon les lois de la physique (Newton), je varie de tant."
• L'Intégration : "Je calcule le point suivant."

La syntaxe

La fonction odeint s'utilise toujours avec 3 arguments principaux.

Syntaxe : Solution = odeint(Modele, Conditions_Initiales, Temps)

• Modele : C'est le nom de la fonction Python qui contient les équations de Newton (PFD). Elle doit renvoyer les dérivées.
• Conditions_Initiales : L'état du système à t=0 (Position et Vitesse de départ).
• Temps : Le tableau NumPy des instants où on veut calculer la position (créé avec linspace).

Le vecteur d'état

Une équation différentielle mécanique est souvent du 2ème ordre (elle contient une accélération x¨). Or, odeint ne sait résoudre que des équations du 1er ordre (y′=f(y)).

L'astuce consiste à grouper la position et la vitesse dans un seul vecteur appelé Vecteur d'État S.

Pour un mouvement 1D (comme la chute libre), on pose : $S = \begin{bmatrix} position \\ vitesse \end{bmatrix}$

On doit alors fournir à Python la dérivée de ce vecteur (S′ ou Sp dans le code) : $S' = \frac{dS}{dt} = \begin{bmatrix} \dot{position} \\ \dot{vitesse} \end{bmatrix} = \begin{bmatrix} vitesse \\ accélération \end{bmatrix}$

Pourquoi cette astuce ?

Cela transforme un problème d'accélération (complexe) en un système simple :

• La dérivée de la position → C'est la vitesse (c'est une définition).
• La dérivée de la vitesse → C'est l'accélération (c'est Newton : ∑F/m).

Cette étape est essentielle pour résoudre des systèmes mécaniques avec odeint.

Résoudre une équation du premier ordre

Prenons l'exemple du boulet de canon ou d'une chute avec frottement.

Supposons l'équation : $\dot{v} = -g - \frac{k}{m}v$.

Voici la structure obligatoire d'un script de résolution :

Passer au second ordre : le pendule

Pour le pendule (et le planeur), on a une position et une vitesse. Le vecteur $S$ a donc deux cases.

L'équation du pendule est : $\ddot{\theta} = -\frac{g}{L}\sin(\theta)$.

Un système complet : le planeur

Pour un planeur, le mouvement se fait en X et en Y. On a donc 4 variables à suivre simultanément : $S = [x, y, v_x, v_y]$

La fonction Sp devra donc retourner 4 valeurs : $S' = [v_x, v_y, a_x, a_y]$

C'est exactement la même logique, mais avec plus de lignes.

Notons sur le Slicing (Rappel)

Quand on récupère la solution S pour le planeur :

• S[:, 0] sera toutes les positions x(t)
• S[:, 1] sera toutes les altitudes y(t)

Pour tracer la trajectoire (la forme du vol), on fera donc :
plt.plot(S[:, 0], S[:, 1]) (x en abscisse, y en ordonnée).

Vérification

Synthèse

• Résoudre une équation différentielle numériquement, c'est avancer pas à pas : on part d'un état connu, on calcule sa dérivée, on en déduit l'état suivant. Aucune formule exacte n'est produite, seulement une suite de points.
• odeint(modele, etat_initial, instants) demande trois choses : une fonction qui rend les dérivées, l'état au départ, et les instants où l'on veut la réponse.
• La fonction modèle ne calcule jamais la solution : elle répond à « dans cet état, à cet instant, comment chaque grandeur varie-t-elle ? ». C'est la traduction directe de l'équation.
• Une équation du second ordre se ramène au premier en empilant position et vitesse dans un vecteur d'état : $S = [x, \dot{x}]$, dont la dérivée est $[\dot{x}, \ddot{x}]$. Le même mécanisme couvre alors les systèmes couplés, en allongeant le vecteur.
• L'ordre des grandeurs dans le vecteur d'état est une convention qu'on choisit, et qu'il faut respecter partout : à la construction, dans le modèle, à la lecture du résultat. La plupart des erreurs de simulation viennent de là.
• Un solveur ne se plaint jamais : il rend des nombres même quand le modèle est faux. L'étalonner d'abord sur un cas dont la solution exacte est connue est la seule façon d'avoir une raison de croire le cas réel.