La méthode d'Euler en Python

8.1 L'idée

On ne sait pas toujours résoudre une équation différentielle exactement. Le programme cite Euler parmi ceux qui ont structuré ces résolutions, aux côtés de Newton, D'Alembert, Lagrange, Cauchy, Clairaut et Riccati. La méthode d'Euler en donne une solution approchée. Elle repose sur la tangente : près d'un point, la courbe est proche de sa tangente. Pour $h$ proche de $0$, $y(x + h) \approx y(x) + h\,y'(x).$ On part du point $(x_0\,;y_0)$ et on avance par petits pas de longueur $h$, chaque fois le long de la droite dont le coefficient directeur est donné par l'équation. On obtient une ligne brisée proche de la courbe de la solution.

8.2 Le cas de $y' = f$

Pour l'équation $y' = f$ avec $y(x_0) = y_0$, on calcule de proche en proche : $x_{k+1} = x_k + h, \qquad y_{k+1} = y_k + h\,f(x_k).$

En langage naturel (« $\leftarrow$ » se lit « prend la valeur »), une étape s'écrit : $y \leftarrow y + h \times f(x)$, puis $x \leftarrow x + h$.

Exemple 40. Pour $y' = \dfrac1x$ avec $y(1) = 0$, la solution exacte est $\ln$. Avec $h = 0{,}1$, le premier pas donne $y_1 = 0 + 0{,}1 \times \dfrac11 = 0{,}1$, le deuxième $y_2 = 0{,}1 + 0{,}1 \times \dfrac{1}{1{,}1} \approx 0{,}190\,9$. Après $10$ pas, on arrive à $x = 2$ avec $y \approx 0{,}718\,8$, alors que $\ln 2 \approx 0{,}693\,1$. Avec $h = 0{,}01$ et $100$ pas, on obtient $0{,}695\,7$ ; avec $h = 0{,}001$ et $1000$ pas, $0{,}693\,4$.

En Python. On range les valeurs dans deux listes, construites par ajouts successifs.

from math import exp

def euler(f, x0, y0, h, n):
    X = [x0]
    Y = [y0]
    for k in range(n):
        Y.append(Y[k] + h * f(X[k]))
        X.append(X[k] + h)
    return X, Y

def f(x):
    return exp(-x**2)

X, Y = euler(f, 0, 0, 0.1, 10)
print(Y[10])

Ce programme traite la fonction $x \mapsto \mathrm{e}^{-x^2}$, qui n'a pas de primitive explicite. Il affiche 0.7778168240731773, valeur approchée de $y(1)$ pour la solution telle que $y(0) = 0$. Avec un pas de 0.01 et $100$ pas, on obtient 0.7499786042621127 ; avec 0.001 et $1000$ pas, 0.747140131778599. La valeur approchée se stabilise vers $0{,}747$.

8.3 Le cas de $y' = ay + b$

La dérivée dépend maintenant de $y$. On remplace $y'$ par $ay + b$ : $x_{k+1} = x_k + h, \qquad y_{k+1} = y_k + h\,(a\,y_k + b).$ La suite $(y_k)$ vérifie $y_{k+1} = (1 + ah)\,y_k + hb$ : c'est une suite arithmético-géométrique.

def euler2(a, b, x0, y0, h, n):
    X = [x0]
    Y = [y0]
    for k in range(n):
        Y.append(Y[k] + h * (a * Y[k] + b))
        X.append(X[k] + h)
    return X, Y

X, Y = euler2(1, 0, 0, 1, 0.01, 100)
print(Y[100])

Ce programme traite $y' = y$ avec $y(0) = 1$, dont la solution est l'exponentielle : il affiche une valeur approchée de $\mathrm{e}$. Tu as vu cette construction en Première. Ici $y_k = (1 + h)^k$.

Ce que l'on observe. Quand on divise le pas par $10$, l'écart est à peu près divisé par $10$. Pour gagner une décimale, il faut dix fois plus de calculs.

Exemple 41 (la tarte). On reprend $T' = -0{,}1\,T + 2$ avec $T(0) = 180$. Avec un pas de $5$ minutes, $y_{k+1} = 0{,}5\,y_k + 10$ et on obtient :

$t$	$0$	$5$	$10$	$15$	$20$
Euler, $h = 5$	$180$	$100$	$60$	$40$	$30$
Solution exacte	$180$	$117{,}04$	$78{,}86$	$55{,}70$	$41{,}65$

L'approximation est grossière. Avec euler2(-0.1, 2, 0, 180, 0.1, 200), on trouve $41{,}44$ à $t = 20$, et avec un pas de $0{,}01$, $41{,}63$ : on se rapproche de la valeur exacte $41{,}65$.

À quoi ça sert, puisqu'on sait résoudre ? Pour les équations du chapitre, on peut comparer à la solution exacte et voir comment se comporte la méthode. Mais elle s'applique aussi aux équations que l'on ne sait pas résoudre : c'est son intérêt principal.