Parlons un peu de la méthode d’Euler. Première question, à quoi ça sert ? C’est une méthode simple qui permet de calculer numériquement la solution d’une équation différentielle du premier ordre. C’est à dire qu’on a une fonction qui exprime la dérivée d’une quantité (comment varie cette quantité dans le temps) et nous, à partir de ça, on va chercher cette quantité.
Petit exemple
Imaginez que vous êtes sur une droite. Vous pouvez juste avancer ou reculer. Comme hier, on se donne une origine sur cette droite pour pouvoir se répérer. Disons que vous êtes à la position $p_0$ à l’instant initial, et que vous ne bougez pas.
Puis, vous bougez à la vitesse $v_1$ jusqu’à l’instant $t_1$, ce qui vous fait alors arriver en $p_1 = p_0 + v_1t_1$. Vous bougez ensuite à la vitesse $v_2$ jusqu’à l’instant $t_2$, ce qui vous donne $p_2 = p_1 + v_2(t_2 - t_1)$. À la vitesse $v_3$ jusqu’à l’instant $t_3$, etc etc. On a donc le schéma suivant :
$$\begin{cases} p(t_0) &= p_0 \newline p(t_{n+1}) &= p(t_n) + v_{n+1}(t_{n+1} - t_n) \end{cases}$$
où $v_{n+1}$ est la vitesse (constante) entre l’instant $t_n$ et $t_{n+1}$. Grosso modo, ça dit que notre position actuelle, c’est la position d’avant à laquelle on ajoute la distance que l’on vient de parcourir.
Voilà. Vous savez quoi ? Si vous avez compris ce qu’on vient de faire là, vous avez compris la méthode d’Euler. C’était pas si compliqué, non ? On va quand même donner la méthode sur le cas général, comme ça après, on pourra l’appliquer à n’importe quoi.
Cas général
Dans l’exemple précédent, nous avions une quantité connue à l’instant initial, la position. Nous savions aussi comment variait cette quantité, grâce à la vitesse. Mathématiquement, nous avions la relation suivante :
$$ p^\prime(t) = v(t) $$
où je note par $^\prime$ la dérivée par rapport au temps. Toujours en gardant cette notation, on prend un cas plus général.
On a une quantité $y$, on connait sa valeur initiale. On sait connait aussi comment elle varie dans le temps (admettons). Disons que cette variation peut s’écrire sous la forme d’une fonction $f(y(t), t)$. On a alors :
$$ y^\prime(t) = f(y(t), t) $$
(pour raccrocher les wagons si besoin, dans l’exemple précédent $y(t)= p(t)$ et $f(y(t), t) = v(t)$)
On se donne un espace temporel, défini comme hier. C’est à dire qu’on a un instant initial $t_0$, un instant final $t_f$ et un pas de temps $dt$. On applique la méthode d’Euler comme on a fait dans l’exemple, ce qui nous donne :
$$\begin{cases} y(t_0) &= y_0 \newline y(t_{n+1}) &= y(t_n) + dt f(y(t_n), t_n) \end{cases} $$
Ainsi, à chaque instant de notre espace temporel, on est capable de calculer la valeur de la quantité $y$.
À un moment, j’ai dit que la vitesse devait être constante dans l’interval $[t_n, t_{n+1}]$. Là aussi, c’est bien évidemment le cas. Pour que ça soit un calcul exact, il faut que la fonction $f(y(t), t)$ soit constante sur chaque interval $[t_n, t_n + dt]$. Sauf qu’en général, c’est pas du tout le cas. Donc la solution à ça (pour le moment), c’est de prendre un pas de temps $dt$ plus petit. Mais après, votre calcul mettra plus de temps… À vous de voir si vous avez besoin d’un calcul précis ou pas.
Un peu de python ?
Pour la dernière parti de ce post, je vous propose de programmer cette méthode. On va faire un truc encore une fois assez simple. On va s’amuser à calculer la vitesse d’un caillou qui tombe.
On se donne donc un petit caillou, d’une masse $m$. On le lève bien haut et on le lâche, sans lui donner de vitesse. On le lâche juste. Donc à l’instant initial, ça vitesse est $v(t = 0) = 0$. Ce petit caillou tout innocent qui est lâché dans le vide n’est soumis qu’à la force de son poids, qui l’entraine vers le bas. Son poids est $P = mg$ où $g$ est l’accélération de la pesenteur. Bien. La seconde loi de Newton nous dit que la somme des forces appliquées à notre petit caillou est égale à la masse du caillou multipliée par son accérélation. Donc on a :
$$ \sum F = P = ma $$
où $a$ est l’accérélation de notre caillou. On notera que l’on peut simplifier par $m$ et que du coup, l’accélération de notre caillou ça devient plus que $g$…
On a dit que nous cherchions la vitesse du caillou. On connait la dérivée par rapport au temps de cette vitesse, qui est l’accélération. On a :
$$ v^\prime(t) = a(t) = g $$
On peut donc appliquer notre méthode d’Euler.
Un peu de code avant pour créer notre espace temporel et les conditions initiales.
>>> from pylab import * >>> g = -9.81 # on définie l’accélération de la pesenteur >>> #on définie notre espace temporel >>> dt = 0.1 >>> tf = 1 >>> times = arange(0, tf, dt) >>> v = 0.0 #la vitesse à l’instant t = 0
Ok. Nous avons défini notre espace temporel. On s’interesse à la vitesse du caillou entre 0 et 1 seconde. Maintenant, on va construire une fonction vprime(v, t), qui va nous donner l’accérélation du caillou. On donne à cette fonction deux paramètres, $v(t_n)$, et $t_n$. Même si on ne se sert pas des deux paramètres dans le calcul, on les donne par convention.
>>> def vprime(v, t): # on se donne une fonction qui définie v’(t) = f(v, t) ... return g
Tout simplement. Et maintenant, le code qui applique la méthode d’Euler.
>>> for t in times: ... print('À l’instant t = {}, la vitesse est v = {}'.format(t, v) ... v = v + dt*vprime(v, t) 'À l’instant t = 0.0, la vitesse est v = 0' 'À l’instant t = 0.1, la vitesse est v = -0.981' 'À l’instant t = 0.2, la vitesse est v = -1.962' 'À l’instant t = 0.3, la vitesse est v = -2.943' 'À l’instant t = 0.4, la vitesse est v = -3.924' 'À l’instant t = 0.5, la vitesse est v = -4.905' 'À l’instant t = 0.6, la vitesse est v = -5.886' 'À l’instant t = 0.7, la vitesse est v = -6.867' 'À l’instant t = 0.8, la vitesse est v = -7.848' 'À l’instant t = 0.9, la vitesse est v = -8.829'
Cool, non ? Oui, c’est cool. Mais, je me dis que ça serait encore plus cool si on connaissait aussi la position de notre caillou… Non ? Comme ça, on pourrait se poser des questions comme « Si je lache mon caillou à 1 mètre du sol, combien de temps lui faut-il pour arriver par terre ? ».
De l’accélération à la position
Dans l’exemple initial, on connaissait la vitesse et on cherchait la position. Là dans l’exemple qu’on vient de faire avec python, on connaissait l’accélération et on cherchait la vitesse. Donc si on met les deux exemples bout à bout, on peut relier l’accélération et la position. Et pour faire ça, on peut le faire avec une méthode imbriquée, telle que je l’ai évoquée ici. On va donc suivre cette méthode et construire un vecteur $Y$ qui contient la vitesse et la position, comme suit :
$$ Y(t) = (v(t), p(t)) $$
On calcule la dérivée de ce vecteur, et on a :
$$ Y^\prime(t) = (a(t), v(t)) = f(Y(t), t)$$
et là, on peut appliquer la méthode d’Euler, ok ? On va appliquer ça tout de suite.
>>> g = -9.81 >>> tf = 1.0 >>> dt = 0.1 >>> times = arange(0, tf, dt) >>> p = 1.0 #On lâche de caillou à 1m de hauteur. >>> v = 0.0 #On ne lui donne pas de vitesse initiale >>> y = array([v, p]) #on construit notre vecteur >>> def yprime(y, t): ... return array([g, y[0]])
On a définit nos conditions initiales, notre fonction qui calcule la dérivée, y’a plus qu’à faire du Euler dessus.
>>> for t in times: ... print('t = {}, v = {}, p = {}'.format(t, y[0], y[1]) ... y = y + dt*yprime(y, t) 't = 0.0, v = 0.0, p = 1.0' 't = 0.1, v = -0.981, p = 1.0' 't = 0.2, v = -1.962, p = 0.9019' 't = 0.3, v = -2.943, p = 0.7057' 't = 0.4, v = -3.924, p = 0.4114' 't = 0.5, v = -4.905, p = 0.019' 't = 0.6, v = -5.886, p = -0.4715' 't = 0.7, v = -6.867, p = -0.0601' 't = 0.8, v = -7.848, p = -0.7468' 't = 0.9, v = -8.829, p = -0.5316'
Et donc là, à vue de nez, on peut voir que si on lâche un caillou à un mètre du sol, l’impact à lieu à environ 0.5 seconde. (Si vous voulez un calcul plus précis, mettez un $dt$ plus petit !)
On s’arrête là pour aujourd’hui. On peut imaginer d’autres trucs cools à faire :
- étendre cet exemple au cas 2D. C’est à dire que le caillou n’est plus laché verticalement, mais jetté en avant, par exemple.
- appliquer la méthode d’Euler sur système masse-ressort.
On a donc vu aujourd’hui la méthode d’Euler. C’était juste pour vous montrer un peu comment on pouvait résoudre des équations différentielles numériquement. Sachez qu’en pratique, dans le monde réel, on n’utilise pas du tout la méthode d’Euler car elle donne très vite des grosses erreurs. Donc demain (ou ce weekend) on parlera de la méthode de Runge-Kutta et on verra que la méthode d’Euler est en fait juste un cas particulier de Runge-Kutta. Et puis, dans un futur un peu plus lointain, on parlera de comment trouver le pas de temps optimal (car pour le moment, faut bien avouer qu’on le met un peu au pif…)
Et voici, comme promis, un peu de code qui résume ce qu’on a fait aujourd’hui. (Le code est un peu plus « compliqué » car c’est le cas 2D… mais ça change pas trop trop non plus…) Pour le fun, j’vous ai aussi mis l’animation graphique ! Enjoy.
Pour vous motiver, un petite animation de ce que vous trouverez dans le code fourni ;) (Ça simule un caillou qu’on jette vers l’avant)
PS : La suite est disponible sur cette page.