SciPy : équations, intégrales et ajustement

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. Ces exemples fonctionnent à l'identique sur une machine.

Une bibliothèque, des sous-modules

SciPy ne s'importe presque jamais en entier. On ouvre le tiroir dont on a besoin :

Tous partagent une habitude qui vaut d'être notée tout de suite : une fonction numérique renvoie un résultat et, souvent, de quoi le juger. Une estimation d'erreur, un indicateur de convergence, une matrice de covariance. Ignorer cette seconde sortie est la faute la plus courante du calcul scientifique.

Trouver un zéro

Beaucoup d'équations d'ingénieur n'ont aucune solution exprimable. Le débit dans une conduite, le point de fonctionnement d'un transistor, l'angle d'équilibre d'un mécanisme : dans tous ces cas, l'équation se met sous la forme $f(x) = 0$ et se résout numériquement.

Le principe le plus répandu est la méthode de Newton : depuis un point, suivre la tangente jusqu'à l'axe des abscisses, et recommencer depuis le point d'arrivée. La tangente est la meilleure droite qui ressemble à la courbe, donc son zéro est une bien meilleure estimation que le point de départ.

En pratique, on n'écrit pas Newton à la main. brentq fait mieux : il combine plusieurs méthodes et garantit de trouver un zéro, à condition de lui donner un intervalle où la fonction change de signe.

La vérification de la dernière ligne n'est pas une coquetterie : elle est le seul élément qui distingue un zéro d'un nombre quelconque. Un solveur numérique rend toujours une valeur, même quand il a convergé vers autre chose que ce qu'on cherchait.

Calculer une intégrale définie

Beaucoup de fonctions n'ont pas de primitive exprimable : $e^{-x^2}$, la plus célèbre, gouverne pourtant toute la statistique. L'aire sous sa courbe se calcule quand même, en la découpant.

Le second résultat est la raison d'être de quad. Une erreur estimée de l'ordre de $10^{-14}$ autorise à faire confiance aux dix premières décimales ; une erreur de l'ordre de $10^{-3}$ signale que la fonction a malmené la méthode, et qu'il faut découper l'intervalle ou en changer.

Ajuster un modèle sur des mesures

Le chapitre sur Matplotlib traçait des relevés et un modèle côte à côte, avec des coefficients choisis à la main. curve_fit les calcule : il cherche les paramètres qui rapprochent le plus le modèle des points, au sens des moindres carrés.

L'argument p0 donne une estimation de départ. Il n'est pas facultatif en pratique : sur un modèle exponentiel, un point de départ absurde fait converger l'algorithme vers un minimum local qui ne ressemble à rien, sans le moindre message d'erreur.

La matrice covariance est la sortie qu'il ne faut pas jeter. La racine de sa diagonale donne l'incertitude sur chaque paramètre, et c'est elle qui permet d'écrire « τ = 0,198 ± 0,004 s » plutôt qu'un nombre nu dont personne ne sait s'il vaut trois décimales ou une.

Prédire une évolution : les équations différentielles

Les lois de Newton, $\sum \vec{F} = m\vec{a}$, relient une accélération aux forces. Ce sont des équations différentielles, et dès que des frottements ou une géométrie un peu riche apparaissent, aucune formule ne donne la trajectoire. Il reste à l'avancer pas à pas.

odeint, du module scipy.integrate, fait exactement cela. Son raisonnement ne contient aucune physique : il ne connaît que la pente.

La fonction modèle ne calcule donc jamais la solution. Elle répond à une seule question, toujours la même : « dans cet état, à cet instant, comment chaque grandeur varie-t-elle ? ». C'est la traduction directe de l'équation, ligne pour ligne.

Syntaxe : solution = odeint(modele, etat_initial, instants)

• modele : la fonction qui renvoie les dérivées, dans l'ordre des grandeurs de l'état.
• etat_initial : la valeur de chaque grandeur à l'instant de départ.
• instants : le tableau des instants où la réponse est souhaitée, fabriqué par linspace.

Le vecteur d'état

Un solveur ne sait résoudre que des équations du premier ordre, de la forme $y' = f(y, t)$. Or la mécanique en produit du second, puisqu'elle porte sur une accélération.

L'astuce consiste à empiler la position et la vitesse dans un seul vecteur, appelé vecteur d'état :

$S = \begin{bmatrix} x \\ \dot{x} \end{bmatrix} \qquad\text{donc}\qquad S' = \begin{bmatrix} \dot{x} \\ \ddot{x} \end{bmatrix} = \begin{bmatrix} \text{vitesse} \\ \text{accélération} \end{bmatrix}$

Les deux lignes de $S'$ n'ont pas le même statut, et c'est ce qui rend l'astuce si simple :

• la dérivée de la position est la vitesse, par définition, donc la première ligne se recopie depuis l'état ;
• la dérivée de la vitesse est l'accélération, que Newton donne : $\sum F / m$.

Une équation du second ordre est ainsi devenue un système de deux équations du premier. Le même mécanisme couvre ensuite n'importe quel système couplé : il suffit d'allonger le vecteur.

Attention : l'ordre des grandeurs dans le vecteur est une convention choisie par celui qui écrit, et elle doit être respectée aux trois endroits : à la construction de l'état initial, dans le modèle, et à la lecture du résultat. La majorité des simulations fausses viennent d'une inversion entre ces trois endroits, et le solveur ne la signale jamais.

Le pas de temps n'est pas celui qu'on croit

Le tableau instants donne les dates où l'on veut lire la réponse, pas celles où le solveur calcule. odeint choisit son propre pas interne, le raccourcit là où la solution varie vite, et l'allonge ailleurs.

La conséquence est utile : demander mille points de sortie au lieu de cent ne rend pas la simulation plus juste, seulement plus détaillée à l'affichage. Ce qui la rend juste est la tolérance du solveur, réglable par rtol et atol. Et quand une solution diverge brutalement alors que le problème est physique, c'est le signe d'un système raide, où une grandeur évolue mille fois plus vite qu'une autre : il faut alors un solveur prévu pour cela, comme solve_ivp(..., method="Radau").

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.

Lire le résultat. La solution rendue par odeint est un tableau à deux dimensions : une ligne par instant, une colonne par grandeur du vecteur d'état. Le découpage NumPy en extrait une colonne entière :

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

La trajectoire, c'est-à-dire la forme du vol, se trace alors en portant l'altitude contre la distance : plt.plot(S[:, 0], S[:, 1]). Porter l'altitude contre le temps donnerait une autre figure, tout aussi juste et répondant à une autre question, celle de la vitesse de descente.