Simulez des systèmes physiques avec la méthode d'Euler

En guise d’introduction, je vous présente l’équation différentielle-type qui nous accompagnera tout au long de ce tutoriel :

$\frac{\mathrm{d}y}{\mathrm{d}t}(t) = f(y(t), t)$

Il s’agit d’une équation différentielle du premier ordre. La variable réelle $t$ représente le temps, compris entre zéro et une valeur limite $T$ ; $y$ est une fonction réelle du temps à déterminer, et $\frac{\mathrm{d}y}{\mathrm{d}t}$ est sa dérivée ; $f$ est une fonction réelle de $y(t)$ et $t$. On note $y_0$ la valeur initiale de $y$, ce qui signifie simplement que $y(0) = y_0$.

Sous des conditions généralement remplies en physique, cette équation différentielle admet une unique solution $y$ pour une condition initiale donnée.¹ Il n’est donc pas vain d’essayer de résoudre cette équation différentielle !

L’équation différentielle ci-dessous est un exemple correspondant à l’équation-type vue précédemment, en prenant $f(y(t), t) = -ay(t) + bt$, où $a$ et $b$ sont des constantes.

$\frac{\mathrm{d}y}{\mathrm{d}t}(t) + ay(t) = bt$

Discrétisation du temps

Il est difficile pour un ordinateur de travailler avec des fonctions continues quelconques. Pour escamoter cet obstacle, il est courant de commencer par discrétiser le temps lors de la résolution numérique des équations différentielles (c’est notamment le cas dans la méthode d’Euler, qui nous intéresse ici). Autrement dit, on ne considère que certains instants particuliers, et on ignore ce qu’il se passe en dehors de ces instants. Formellement, on forme une suite croissante d’instants $t_k$ pris dans l’intervalle $[0, T]$.

$(t_k)_{k \in \mathbb{N}}$

Ces instants peuvent a priori être choisis librement (par exemple toutes les secondes), mais la méthode de résolution numérique utilisée oriente en général ce choix. Notez que j’ai indexé cette suite dans $\mathbb{N}$, bien que les instants soient en nombre fini, car on peut en prendre autant que l’on souhaite.

Si l’on connaît une solution exacte $y$ de l’équation différentielle, elle prend évidemment la valeur $y(t_k)$ à l’instant $t_k$. La résolution numérique consiste à calculer approximativement les $y(t_k)$, et obtenir ainsi une solution numérique, qu’on espère fidèle à la solution exacte. La manière précise de calculer les $y_k$ dépend de la méthode de résolution utilisée. Formellement, on construit une suite de valeurs $y_k$ approchant les $y(t_k)$.

$\forall k \in \mathbb{N},~ y_k \approx y(t_k)$

Pour la représentation graphique, on considère que la solution approchée est linéaire entre deux instants consécutifs. Notez bien que cette hypothèse n’est utilisée que pour tracer graphiquement la fonction, et n’est en aucun cas une hypothèse nécessaire pour utiliser la méthode d’Euler.

Approximation d’Euler

En choisissant de discrétiser le temps, nous sommes passés du monde des fonctions du temps à celui des suites numériques. Derrière cette apparente simplification se cache un obstacle à surmonter. En effet, le problème à résoudre est une équation différentielle, qui fait intervenir la notion de dérivée. Or cette notion fait intervenir la notion de limite (la dérivée est la limite, si elle existe, du taux d’accroissement), qu’on ne peut pas définir pour les suites.

Pour lever le problème, il faut adapter intelligemment l’opérateur de dérivation aux suites numériques ; c’est là que la méthode d’Euler entre véritablement en jeu. Dans sa variante dite explicite, que nous traitons ici, elle consiste à dire que le taux d’accroissement entre un instant $t_k$ et un instant ultérieur $t_{k+1}$ est égal à la dérivée à l’instant $t_k$. Autrement dit, on approxime le taux d’accroissement à l’instant $t_{k}$ par :

$\frac{y(t_{k+1}) - y(t_k)}{t_{k+1} - t_k} \approx \frac{\mathrm{d}y}{\mathrm{d}t}(t_k)$

La figure ci-dessous illustre cette approximation.

On peut déjà avoir l’intuition que l’écart entre deux instants ne saurait être trop important, au risque d’avoir une dérivée trop éloignée de la réalité ! Essayez de visualiser cela sur la figure précédente : en rapprochant $t_{k+1}$ de $t_k$, l’accroissement approximé se rapproche de l’accroissement exact à mesure que la corde se rapproche de la tangente.

Discrétisation de l’équation différentielle

Nous avons maintenant tout le nécessaire pour transformer notre équation différentielle en suite récurrente, et compléter ainsi le passage du monde continu au monde des suites numériques.

Je vous rappelle ci-dessous l’équation différentielle que l’on cherche à résoudre numériquement sur l’intervalle $[0, T]$ :

$\frac{\mathrm{d}y}{\mathrm{d}t}(t) = f(y(t), t)$

En particulier, cette équation différentielle est vraie pour $t = t_k$, de sorte que si $y$ est solution, $y(t_k)$ vérifie l’équation suivante :

$\frac{\mathrm{d}y}{\mathrm{d}t}(t_k) = f(y(t_k), t_k)$

On applique l’approximation d’Euler pour faire disparaître la dérivée et obtenir une nouvelle équation, qui n’est vérifiée qu’approximativement par la solution exacte :

$\frac{y(t_{k+1}) - y(t_k)}{t_{k+1} - t_k} \approx f(y(t_k), t_k)$

Bien qu’elle ne soit qu’approximative, cette équation est très utile, car elle donne une recette de cuisine pour fabriquer une solution approchée. En effet, on peut décider que la solution approchée vérifie cette égalité :

$\frac{y_{k+1} - y_k}{t_{k+1} - t_k} = f(y_k, t_k)$

Ceci fait, on en déduit une relation de récurrence en gardant seulement $y_{k+1}$ à gauche et en transférant le reste dans le membre de droite. On appelle cette relation le schéma numérique d’Euler explicite :

$y_{k+1}= y_k + (t_{k+1} - t_k)f(y_k, t_k)$

$y_k$ est issu de l’étape précédente, les $t_k$ sont choisis librement, et $f$ est une fonction connue. On peut ainsi calculer $y_{k+1}$ de proche en proche en partant de la valeur initiale $y_0$. Nous avons transformé une équation différentielle en une simple suite définie par récurrence ! Il se trouve que c’est une formulation très avantageuse pour un calcul par ordinateur, car ils sont justement conçus pour effectuer des opérations successives. C’est ce que nous ferons en passant à la pratique dans la partie suivante !

Position du problème

Considérons la chute verticale vers le bas d’un objet dans l’air, comme un caillou par exemple. L’objet est attiré par gravité vers le sol, mais il est aussi freiné par l’air, et il sera d’autant plus freiné que sa vitesse sera grande. On peut donc intuitivement penser que sa vitesse augmentera progressivement, avant d’atteindre une vitesse limite.

Pour vérifier cela, modélisons notre problème. L’objet a une masse $m$, et est suffisamment dense pour qu’on néglige la poussée d’Archimède. Il est soumis au champ de pesanteur d’intensité $g$ ; il subit également une force de frottement exercée par l’air proportionnelle à la vitesse d’un facteur $f$ et opposée au mouvement. On considérera sa vitesse initiale nulle. En appliquant le principe fondamental de la dynamique à l’objet, on obtient une équation différentielle régissant la vitesse $v$ de l’objet.

$\frac{\mathrm{d}v}{\mathrm{d}t}(t) + \frac{f}{m}v(t) = g$

Remarques préliminaires

Sans résoudre analytiquement l’équation différentielle, il est possible de confirmer notre intuition concernant la vitesse limite. En effet, si les frottements équilibrent la gravité, on peut écrire la relation suivante concernant la vitesse d’équilibre $v_{eq}$ :

$\frac{f}{m}v_{eq} = g$

Cette relation a pour conséquence d’annuler la dérivée :

$\frac{\mathrm{d}v_{eq}}{\mathrm{d}t} = g - \frac{f}{m}v_{eq} = 0$

Une dérivée nulle signifie qu’il n’y a plus d’évolution : la vitesse d’équilibre entre frottements et gravité est aussi la vitesse atteinte à l’état final.

Ce petite calcul montre ainsi qu’en plus d’avoir des informations sur le point de départ (la vitesse initiale ainsi que la pente à l’origine dans notre cas), on en a aussi sur le point d’arrivée (la vitesse finale). Il nous manque juste la transition entre ces deux moments, que l’on peut obtenir approximativement en résolvant l’équation numériquement avec la méthode d’Euler explicite.

Équation de récurrence

On applique la méthode d’Euler pour résoudre l’équation différentielle. Concrètement, cela revient à dire que la solution approchée, la solution numérique, vérifie l’équation obtenue en appliquant l’approximation d’Euler.

$\frac{v_{k+1} - v_k}{t_{k+1} - t_k} + \frac{f}{m}v_k = g$

On isole $v_{k+1}$ à gauche et on transfère le reste dans le membre de droite pour obtenir une relation de récurrence :

$v_{k+1} = v_k + (t_{k+1} - t_k)\left(g - \frac{f}{m} v_k \right)$

On reconnaît une forme similaire à celle de la section précédente, où l’on aurait $f(t,v) = g - \frac{f}{m} v$.

Résolution numérique

On utilise une suite d’instants $t_k$ successifs pour la résolution numérique, et il est courant d’utiliser des pas de temps constants. Cela consiste à choisir les $t_k$ tels que leur écart soit constant et égal à $\Delta t$.

$\forall k,~t_{k+1} - t_k = \Delta t$

On introduit cette notation dans l’équation.

$v_{k+1} = v_k + \Delta t \left( g - \frac{f}{m} v_k \right)$

On prendra les valeurs numériques :

• $g = 10~\mathrm{m.s^{-2}}$,
• $\frac{f}{m} = 1~\mathrm{s^{-1}}$,
• et $\Delta t = 500~\mathrm{ms}$. Cette valeur est acceptable, mais n’offre pas une grande précision, comme nous allons le voir plus loin.

Nous sommes bons pour calculer ! Les ordinateurs font ça très bien, et j’ai écrit à cette fin un court programme Python, qui se trouve dans le bloc masqué ci-dessous. Si vous savez programmer, je vous conseille d’essayer d’écrire votre propre programme avant toute chose, afin de mettre un peu les mains dans le cambouis. Si vous ne savez pas programmer, il est possible d’obtenir un résultat similaire avec un tableur comme Excel. Dans les deux cas, il est intéressant d’avoir un outil pour tester facilement différents paramètres, et ce le sera encore plus dans la partie suivante.

"""
Simulation numérique avec la méthode d'Euler explicite
"""

# Paramètres
t_0 = 0      # Instant initial (s)
y_0 = 0      # Valeur initiale (m/s)
dt = 500e-3  # Pas de temps (s)
t_fin = 10   # Instant final (s)

# Fonction f associée à l'équation différentielle
def f(y, t):
    g = 10  # Accélération de la pesanteur (m/s^2)
    k_sur_m = 1 # f/m (Hz)
    f = g - k_sur_m * y
    return f

# Initialisation de la simulation
t_k = t_0
y_k = y_0

# Initialisation des listes pour le stockage du résultat
t = [t_k]
y = [y_k]

# Déroulement de la simulation
while t_k < t_fin: 
    # Calcul de la solution à l'instant y[k+1]
    y_k = y_k + dt * f(y_k, t_k)
    # Avancement du temps vers l'instant t[k+1]
    t_k = t_k + dt
    # Stockage du résultat
    t.append(t_k)
    y.append(y_k)

# Affichage
print(t)
print(y)

Je consigne les résultats que m’a donnés le programme dans le tableau suivant. Il est redondant afin de bien montrer que le $v_{k+1}$ d’une étape devient le $v_k$ de l’étape suivante.

Comme on peut s’y attendre, la vitesse atteint une valeur limite de 10 m/s après avoir augmenté progressivement¹… Mais à quel point la solution numérique est-elle proche de la réalité ?

Comparaison avec la solution exacte

Il se trouve (comme c’est commode !) que notre équation différentielle admet une solution exacte, qui est la suivante :

$v(t) = \frac{gm}{f} + \left(v(0) - \frac{gm}{f}\right)e^{-\frac{f}{m}t}$

On peut voir clairement la solution exacte et la solution numérique sur le graphique suivant.

Et le verdict est sans appel : la solution numérique ne colle pas à la solution théorique ! Cela est dû au grand pas de temps que nous avons choisi. Comment alors obtenir une bonne précision ? Réduire le pas est probablement une solution… Nous traitons cette question plus en avant dans l’extrait suivant.

Précision en fonction du pas de temps

Puisque la seule chose imprécise est l’approximation de la dérivée, il faut l’améliorer pour gagner en précision. Sans changer de formule d’approximation, il n’y a pas beaucoup de possibilités : il faut réduire le pas de temps. L’idée est de se rapprocher de la définition formelle de la dérivée, qui est la limite (si elle existe) du taux d’accroissement quand le pas de temps tend vers zéro.

$\lim_{\Delta t \rightarrow 0} \frac{y(t+\Delta t) - y(t)}{\Delta t} = \frac{\mathrm{d}y}{\mathrm{d}t} (t)$

La figure suivante illustre cela sur l’exemple de l’extrait précédent. La réduction du pas de temps permet de s’approcher plus près de la solution exacte, jusqu’à la toucher, en théorie, pour un pas de temps nul. On dira alors que la méthode est convergente. On peut aussi s’intéresser à la vitesse de convergence en fonction du pas de temps : je fournis un lien traitant de ce sujet dans la conclusion.

Stabilité en fonction du pas de temps

La méthode d’Euler n’est pas infiniment précise : elle introduit des erreurs d’approximation. Mais une fois qu’une erreur s’est glissée dans le processus, ce qui arrive dès la première étape de calcul, que devient-elle ? Elle peut être soit amplifiée, on parlera de méthode instable, soit au contraire atténuée, et on dira que la méthode est stable.

Dans l’exemple précédent, la solution exacte et la solution numérique se rapprochaient jusqu’à n’être plus discernables à l’œil nu, ce qui signifie que les erreurs s’atténuaient et qu’on était dans un cas stable. Pourtant, la méthode d’Euler n’est pas toujours stable, comme nous allons le voir en analysant plus finement l’exemple.

On note $v(t_k)$ les valeurs prises par la solution exacte de l’équation différentielle. En effectuant un développement limité en $t_k$¹, on obtient : $v(t_{k+1}) = v(t_k+\Delta t) = v(t_k) + \Delta t \frac{\mathrm{d}v}{\mathrm{d}t}(t_k) + e_k$

où $e_k$ est le terme d’erreur du développement limité² ; c’est-à-dire la différence entre la valeur exacte de $v(t_{k+1})$ et la somme des deux premiers termes. Ce terme vérifie $e_k = o(\Delta t)$, c’est-à-dire qu’il est négligeable devant le pas de temps quand celui-ci tend vers zéro. On peut ensuite éliminer la dérivée grâce à l’équation différentielle :

$v(t_{k+1}) = v(t_{k}) + \Delta t \left(g - \frac{f}{m}v(t_{k})\right) + e_k$

Puisqu’à chaque étape, le calcul se base sur l’étape précédente, les termes d’erreurs du développement limité s’accumulent, et il en résulte une erreur globale d’approximation entre la solution approchée et la solution exacte. À une étape $k$ du calcul, on la notera $\delta_k$. Nous avons ainsi $v_k = v(t_k) + \delta_k$. À l’étape suivante, l’erreur se propage, et nous avons donc similairement $v_{k+1} = v(t_{k+1}) + \delta_{k+1}$. Dans le schéma d’Euler, on effectue alors les substitutions suivantes :

• on remplace $v_{k}$ par $(v(t_k) + \delta_k)$,
• et on remplace $v_{k+1}$ par $(v(t_{k+1}) + \delta_{k+1})$.

On obtient alors la relation suivante :

$v(t_{k+1}) + \delta_{k+1} = v(t_{k}) + \delta_k + \Delta t \left(g -\frac{f}{m} (v(t_{k}) + \delta_k)\right)$

Pour se fixer un peu mieux les idées, les différentes erreurs sont illustrées sur la figure ci-dessous.

En soustrayant les deux équations précédentes, on obtient l’équation régissant l’erreur d’approximation. $\delta_{k+1} = \delta_k \left( 1 - \frac{f}{m} \Delta t \right) - e_k$

Par simplification, on laisse tomber $e_k$, tout en conservant l’égalité³. $\delta_{k+1} = \delta_k \left(1 - \frac{f}{m} \Delta t \right)$

Il s’agit d’une récurrence géométrique, qu’il est facile d’étudier. Quelles sont les conditions à remplir pour que l’erreur reste bornée et donc que la méthode soit stable ? Il faut que la suite géométrique soit convergente, ce qui se traduit par une condition sur sa raison : $\left|1 - \frac{f}{m} \Delta t \right| < 1$

Les variables étant positives, la valeur absolue ne peut pas devenir supérieure à 1, et on ne s’intéresse qu’au cas négatif. Il faut dans ce cas : $1 - \frac{f}{m} \Delta t > -1$

En transformant l’inéquation, cela se traduit en condition sur $\Delta t$ ; il nous faut pour assurer la stabilité : $\Delta t < 2\frac{m}{f} = 2~\mathrm{s}$

La figure ci-dessous illustre le comportement pour différentes valeurs de $\Delta t$.

• Pour un pas de 0,6 seconde, la méthode est stable.
• Pour un pas de 1,5 seconde, la méthode est stable mais présente des oscillations numériques, car $1 - \frac{f}{m} \Delta t$ est négatif. Ce phénomène est embêtant, car les simulations ne reflètent même plus la réalité de manière qualitative !
• Pour un pas de 2,0 secondes, on se trouve à la limite de stabilité, la solution oscille sans jamais n’être ni amplifiée, ni attenuée.
• Pour un pas de 2,5 secondes, la solution diverge totalement.