Tous droits réservés

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

Initiation à la simulation numérique

Publié :

En sciences, de nombreux phénomènes sont décrits par des systèmes d'équations différentielles, mais dans certains cas, il est impossible ou peu intéressant de les résoudre de manière exacte. Il s'agit alors de laisser de côté les méthodes formelles pour travailler directement avec des nombres : c'est que qu'on appelle la résolution numérique d'équations différentielles. Résoudre numériquement des équations permet de simuler des systèmes physiques, ce qui en fait un outil puissant pour les sciences et l'ingénierie.

L'objectif de ce tutoriel est de vous initier à la simulation numérique à travers la méthode la plus élémentaire qu'il soit : la méthode d'Euler. Il est destiné aux gens ayant déjà été introduits aux équations différentielles de manière formelle et curieux de voir comment les résoudre autrement.

Pour suivre au mieux ce tutoriel,

  • il est essentiel de savoir ce qu'est une dérivée, une équation différentielle, et une suite numérique ;
  • il est recommandé de savoir résoudre des équations différentielles basiques, et d'avoir quelques notions de physique pour mieux apprécier les exemples.

Présentation de la méthode d'Euler

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

$$ \left.\frac{\mathrm{d}y}{\mathrm{d}t}\right|_t = f(t,y(t)) $$

Il s'agit d'une équation différentielle du premier ordre. $t$ est une variable réelle qui représente le temps, comprise entre zéro et une valeur limite $T$ ; $y$ est une fonction réelle du temps à déterminer ; et $f$ est une fonction réelle de $t$ et de $y(t)$. On note $y_0$ la valeur initiale de $y$.

$$ y(0) = y_0 $$

La notation $\left.\frac{\mathrm{d}y}{\mathrm{d}t}\right|_t$ est synonyme de $\frac{\mathrm{d}y}{\mathrm{d}t}(t)$. Autrement dit, c'est la valeur de la dérivée par rapport au temps en $t$.

Résoudre notre équation différentielle consiste donc à déterminer les valeurs prises sur un certain intervalle par une certaine fonction présentant des contraintes sur son évolution. Sous des conditions généralement remplies en physique, cette équation différentielle admet une unique solution $y$ pour une condition initiale donnée.1 Il n'est donc pas vain d'essayer de résoudre cette équation différentielle !

Un exemple d'équation différentielle suivant le modèle précédent est donné ci-dessous, avec $f(t,y(t)) = -y(t) + at$, où $a$ est une constante réelle dont la dimension physique est telle que $y(t)$ et $at$ soient de même dimension.

$$ \left.\frac{\mathrm{d}y}{\mathrm{d}t}\right|_t + y(t) = at $$

Faisons-nous discrets

Les ordinateurs ne peuvent manipuler des entités qu'en nombre fini. Ainsi, la première chose à faire lorsqu'on veut résoudre numériquement une équation différentielle, c'est de discrétiser le temps : au lieu d'avoir une variable continue, on ne considère que des instants particuliers. Plus formellement, on forme une suite d'instants $t_k$ pris dans l'intervalle $[0, T]$.

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

On ne considère que certains instants sur la flèche du temps.

Ces instants sont choisissables librement. Par exemple, la suite pourrait être 0s, 1s, 2s, 3s, etc. Dans une utilisation pratique, on ne considère qu'un nombre fini d'instants, mais j'ai tout de même indexé cette suite dans $\mathbb{N}$, car on peut théoriquement en prendre autant que l'on souhaite.

Aux instants $t_k$, la solution exacte prend naturellement la valeur $y(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. Plus formellement, on construit une suite de valeurs approchées des $y(t_k)$, qu'on notera $y_k$, et qui dépend de la méthode employée.

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

On cherche des valeurs approchées $y_k$ de $y(t_k)$ à des instants bien définis.

Nous nous sommes maintenant débarrassés de la continuité de nos variables ! Malheureusement, la dérivée est définie comme la limite (si elle existe) du taux d'accroissement, et une telle limite n'a pas de sens pour un temps discret. Il est ainsi temps d'adapter l'opérateur de dérivation afin de pouvoir traduire l'équation différentielle en une relation de récurrence ; c'est là que la méthode d'Euler entre en jeu. Ce tutoriel traite de la variante qu'on appelle méthode d'Euler explicite.

La méthode consiste à approximer la dérivée par le taux d'accroissement entre l'instant actuel et l'instant suivant. On approxime ainsi la dérivée à l'instant $t_{k}$ par :

$$\left. \frac{\mathrm{d}y}{\mathrm{d}t}\right|_{t = t_k} \approx \frac{y(t_{k+1}) - y(t_k)}{t_{k+1} - t_k} $$

En pratique, on ne dispose que des valeurs approchées, qu'on substitue donc aux valeurs exactes dans la formule.

$$ \left.\frac{\mathrm{d}y}{\mathrm{d}t}\right|_{t = t_k} \approx \frac{y_{k+1} - y_k}{t_{k+1} - t_k} $$

On peut améliorer l'approximation en rapprochant $t_{k+1}$.

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'approximation se rapproche de la véritable tangente.

Discrétisation de l'équation différentielle

En appliquant l'approximation précédente à l'équation différentielle-type, on obtient, à l'instant $t_k$ :

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

En gardant seulement $y_{k+1}$ à gauche et en transférant le reste dans le membre de droite, on obtient une relation de récurence qui est le schéma numérique d'Euler.

$$ y_{k+1} = y_k + (t_{k+1} - t_k)f(t_k,y_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 donc 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 ces machines sont justement conçues pour effectuer des opérations successives. C'est ce que nous ferons en passant à la pratique dans la partie suivante !


  1. Pour plus de détails, intéressez-vous au théorème de Cauchy-Lipschitz

Exemple travaillé : chute avec frottements fluides

Position du problème

Considérons la chute verticale vers le bas d'un objet dans l'air (une crotte de pigeon par exemple1). L'objet est attiré par gravité vers le sol, mais il est aussi freiné par l'air. Il sera d'autant plus freiné que sa vitesse sera grande. On peut donc intuitivement penser que sa vitesse va progressivement augmenter, puis 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.

$$ \left.\frac{\mathrm{d}v}{\mathrm{d}t}\right|_{t} + \frac{f}{m}v(t) = g $$

On choisit de résoudre cette é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. On remplace ainsi la dérivée par son approximation.

$$ \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)(g - \frac{f}{m} v_k) $$

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(g - \frac{f}{m} v_k) $$

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 pour cela 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.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
"""
Simulation numérique avec la méthode d'Euler explicite
"""


# Paramètres

# Conditions initiales
t_0 = 0
y_0 = 0

# Fonction f correspondant à l'équation différentielle
def f(t_k,y_k):
    g = 10
    k_sur_m = 1
    f_k = g - k_sur_m*y_k
    return f_k

dt = 500e-3 # Pas de temps
t_fin = 5 # Instant final


# Simulation

# Intialisation
t_k = t_0
y_k = y_0

# Stockage du résultat dans des listes
t = [t_k] 
y = [y_k] 

# Déroulement de la simulation
while t_k < t_fin: 
    # Avancement d'un pas de temps
    y_k = y_k + dt*f(t_k,y_k)
    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.

k $t_k$ $v_k$ $v_{k+1}$
0 0,0 s 0 m/s 5,0 m/s
1 0,5 s 5,0 m/s 7,5 m/s
2 1,0s 7,5 m/s 8,75 m/s
3 1,5 s 8,75 m/s 9,38 m/s
4 2,0 s 9,38 m/s 9,69 m/s
5 2,5 s 9,69 m/s 9,84 m/s
6 3,0 s 9,84 m/s 9,92 m/s
7 3,5 s 9,92 m/s 9,96 m/s
8 4,0 s 9,96 m/s 9,98 m/s
9 4,5 s 9,98 m/s 9,99 m/s

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

Comparaison avec la solution exacte

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

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

On peut voir clairement la solution exacte et la solution approchée sur le graphique suivant.

La simulation ne correspond que vaguement à la théorie.

Et le verdict est sans appel : la solution approchée ne colle pas à la solution théorique ! Cela est dû au grand pas de résolution 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.


  1. Un physicien de passage à Paris trouverait ça tout à fait fascinant. 

  2. La limite exacte est effectivement 10 m/s. On peut le démontrer à la physicienne en considérant une solution stationnaire ($\left.\frac{\mathrm{d}v}{\mathrm{d}t}\right|_{t=\infty} = 0 $) à l'infini, ou rigoureusement en résolvant formellement l'équation différentielle. Vous pouvez vérifier cela à titre d'exercice ! 

Précision, convergence et stabilité

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.

$$ \left. \frac{\mathrm{d}y}{\mathrm{d}t} \right|_t = \lim_{\Delta t \rightarrow 0} \frac{y(t+\Delta t) - y(t)}{\Delta 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 : je fournis un lien traitant de ce sujet dans la conclusion.

La précision diminue avec l'augmentation du pas de temps.

En théorie, faire tendre le pas de temps vers zéro permet d'avoir une précision infinie. En pratique, les nombres à virgule flottante (une approximation des nombres réels) utilisés par les ordinateurs ont une précision finie, et lorsque le pas de temps est minuscule, les erreurs de calcul font perdre en précision. Il y a donc un pas de temps minimal dépendant de la précision des nombres sur la machine effectuant les calculs.

Stabilité en fonction du pas de temps

Nous avons vu que la méthode d'Euler n'est pas infiniment précise : elle introduit des erreurs d'approximation, qui font que l'on ne tombe pas parfaitement sur la solution exacte. 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, on dira que la méthode est stable.

Dans l'exemple précédent, les erreurs s'atténuaient, puisque la solution exacte et la solution numérique se rapprochaient jusqu'à n'être plus discernables à l’œil nu. Pourtant, la méthode d'Euler n'est pas inconditionnellement 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$1, on obtient :

$$ v(t_{k+1}) = v(t_k+\Delta t) = v(t_k) + \Delta t \left.\frac{\mathrm{d}v}{\mathrm{d}t}\right|_{t=t_k} + e_k $$

$e_k$ est le terme d'erreur du développement limité2 ; 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 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 (g - \frac{f}{m}v(t_{k})) + 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 remplace donc $v_{k}$ par $(v(t_k) + \delta_k)$ et $v_{k+1}$ par $ (v(t_{k+1}) + \delta_{k+1}) $ :

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

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

$e_k$ est l'erreur faite sur une étape en partant d'une valeur exacte, $\delta_k$ est l'erreur accumulée.

En soustrayant les deux équations précédentes, on obtient l'équation régissant l'erreur d'approximation.

$$ \delta_{k+1} = \delta_k (1 - \frac{f}{m} \Delta t) - e_k $$

Par simplification, on laisse tomber $e_k$, tout en conservant l'égalité3.

$$ \delta_{k+1} = \delta_k (1 - \frac{f}{m} \Delta t) $$

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.

La méthode d'Euler devient instable pour des pas de temps trop grands.


  1. Mathématiquement, il suffit que $v$ soit dérivable deux fois en $t_k$, et d'appliquer la formule de Taylor. Dans le modèle pris pour l'exemple, c'est le cas. 

  2. C'est aussi ici l'erreur de consistance, c'est-à-dire l'erreur que la méthode d'Euler introduit à l'étape suivante en partant d'une valeur juste. C'est un outil théorique, puisque cette erreur s'introduit dès la première étape de calcul, ce qui fait qu'aucune étape ultérieure ne se fonde sur une valeur exacte. Pire encore, la première valeur est, en général, issue du monde réel et donc déjà entachée d'erreurs de mesure ! 

  3. Ce n'est pas du tout rigoureux ! On peut cependant raisonnablement supposer que les $e_k$ sont bornés. Dans ce cas-là, si la suite $(\delta_k)$ ne diverge pas, elle sera aussi bornée. 


Vous avez découvert dans ce tutoriel la résolution numérique d'équations différentielles avec la méthode d'Euler explicite.

Sur un exemple, nous avons vu comment augmenter la précision de la méthode en réduisant précautionneusement le pas de temps, et comment un pas de temps trop grand pouvait mener à une instabilité de la résolution numérique. On peut ainsi obtenir un résultat satisfaisant, mais cela n'empêche pas de faire mieux avec d'autres méthodes !

Si vous savez programmer, un bon exercice d'appropriation est d'adapter le petit programme que j'ai donné pour résoudre une autre équation différentielle, et d'observer le comportement de la solution numérique lorsqu'on fait varier les différents paramètres.

Si vous voulez aller voir plus loin, je vous conseille d'abord de vous intéresser à un traitement plus mathématique de la méthode d'Euler, qui formalise les intuitions que nous avons eues sur l'exemple. Mais si vous n'aimez pas les mathématiques, passez votre chemin.

Ensuite, rien de tel que d'enchaîner sur la méthode d'Euler implicite (aussi simple, mais plus stable), ainsi que les méthodes de Runge-Kutta (avec une approximation plus précise de la dérivée, à pas de temps égal), puis éventuellement sur les méthodes à pas variable (on avance vite quand le signal varie peu et lentement quand il varie rapidement).

Illustration du tutoriel : détail d'un portrait d'Euler par Emanuel Handmann, 1753.

Je remercie vivement Vayel, Gabbro, Algue-Rythme, adri1 et Holosmos pour leurs relectures et conseils.

2 commentaires

tuto tres sympas :) on l'avait vu en cours et c'est super pratique pour resoudre toute sorte d'équation différentielle ( un petit script python, matplotlib et ca roule tout seul ). On a vu une version améliorer aussi, avec une méthode de prédicateur/correcteur : pour l'équation différentielle

$$\dfrac{dx}{dt}=f(x(t),t)$$
on fait une premiere approximation a partir de la dérivée, sans partir pour autant vers Runge-Kutta:
$$ x_{pr}(t_{k+1}) = x(t_k) + \Delta t *f(x(t_k),t_k) $$
et on fait la moyenne avec la valeur classique:
$$ x(t_{k+1}) = x(t_k) + \frac 12 (f(x(t_k),t_k)+f(x_{pr}(t_{k+1}),t_{k+1}) $$

De mémoire on a utilisé cette méthode sur la trajectoire d'un électron dans un champ $\overrightarrow{B}$; La méthode classique donnait une spirale, celle avec prédiction un cercle qui était la trajectoire réelle

Édité par hobi1

+0 -0
Vous devez être connecté pour pouvoir poster un message.
Connexion

Pas encore inscrit ?

Créez un compte en une minute pour profiter pleinement de toutes les fonctionnalités de Zeste de Savoir. Ici, tout est gratuit et sans publicité.
Créer un compte