euler — Exemples d’application de la méthode d’Euler

La méthode d’Euler est une méthode permettant de calculer des solutions approchées d’équations différentielles.

Cet outil fournit différents exemples de mise en œuvre de cette méthode.

Usage

Le but de cet outil est de montrer la simplicité de la mise en œuvre de cette méthode. Pour cette raison, les programmes sont réduits à leur strict minimum, et ne contiennent pas d’options. Pour les appeler, lancer l’une des commandes suivantes

euler chute
euler proiepredateur
euler ressort
euler satellite

Principe et Algorithme

La méthode d’Euler est utilisée pour trouver une solution approchée d’une équation différentielle. Dans toute cette page, elle permet de tracer point par point une solution approchée de la solution.

Problème

Considérons l’équation différentielle suivante :

\[\begin{split}\left\{\begin{array}{rcl} f(0) &=& 1\\ f'(x) &=& 0,5\times f(x)\\ \end{array}\right.\end{split}\]

Pour quelqu’un ayant étudié un tout petit peu les équations différentielles, trouver la solution exacte est simple. Mais il est fréquent, en modélisant des problèmes réels, de devoir résoudre des équations différentielles complexes, dont on ne connaît pas de solution exacte.

Supposons que nous ne connaissions pas la solution exacte. Comment tracer la courbe d’une solution approchée ?

Approximation

Prenons la courbe représentative d’une fonction, et sa tangente en une absicsse \(a\) donnée, comme sur l’exemple ci-dessous.

System Message: WARNING/2 ( \begin{tikzpicture}[ultra thick, spy using outlines={circle,lens={scale=3}, size=2cm, connect spies} ] \begin{scope} \clip (-.5, -.5) rectangle (2.5, 3.5); \draw[-latex] (-.5, 0) -- (2.5, 0); \draw[-latex] (0, -1) -- (0, 3); \draw[color=gray, dotted, step=.5](-.5, -1) grid (2.5, 3); \draw (0, 0) node[below left]{$\mathcal{O}$}; \draw[smooth,domain=-0.5:2.5] plot(\x,{\x*(\x-1)*(\x-2)+1}); \draw[blue,smooth,domain=1.6:2.5] plot(\x, {2*\x-3}); \draw[dashed] (2, 0) node[below]{$a$} -- (2, 1); \end{scope} \spy[red] on (2, 1) in node [fill=white] at (1.2, 2.7); \end{tikzpicture} )

! Package pgfkeys Error: I do not know the key '/tikz/spy using outlines', to w\nhich you passed 'circle,lens={scale=3}, size=2cm, connect spies', and I am goin\ng to ignore it. Perhaps you misspelled it.\n\nSee the pgfkeys package documentation for explanation.\nType H <return> for immediate help.\n ... \n \nl.16 ...e,lens={scale=3}, size=2cm, connect spies} ]\n \n! Undefined control sequence.\nl.27 \\spy\n [red] on (2, 1) in node [fill=white] at (1.2, 2.7);\n[1{/var/lib/texmf/fonts/map/pdftex/updmap/pdftex.map}]\n(./tikz-7cf92799f3a6ce8d4203f8656645c537b2ed2d66.aux))\n(see the transcript file for additional information)\n 425 words of node memory still in use:\n 5 hlist, 2 vlist, 2 rule, 2 glue, 1 kern, 5 attribute, 48 glue_spec, 5 attri\nbute_list, 1 write nodes\n avail lists: 2:529,3:308,4:17,5:24,6:89,7:289,8:2,9:42\n</usr/share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmmi12.pfb></usr/\nshare/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmsy10.pfb>\nOutput written on tikz-7cf92799f3a6ce8d4203f8656645c537b2ed2d66.pdf (1 page, 18\n310 bytes).\nTranscript written on tikz-7cf92799f3a6ce8d4203f8656645c537b2ed2d66.log.\n

On remarque qu’au voisinnage de l’abscisse \(a\), la corube et sa tangente sont presque confondues. C’est l’approximation faite par la méthode d’Euler.

https://upload.wikimedia.org/wikipedia/commons/9/9e/Integration_x_div_2.png

Fig. 8 Intégration de la fonction x/2 par la méthode d’Euler, Pdebart (source).

Premier exemple

Reprenons l’équation différentielle précédente.

\[\begin{split}\left\{\begin{array}{rcl} f(0) &=& 1\\ f'(x) &=& 0,5\times f(x) \\ \end{array}\right.\end{split}\]

La solution est la fonction \(x\mapsto e^{0,5\times x}\). Supposons que nous ne connaissons pas cette réponse, et résolvons cette équation sur l’intervalle \([0;1]\) avec la méthode d’Euler, en traçant la courbe sur cet intervalle, avec un pas \(k=0,1\).

  • Tout d’abord, \(f(0)=1\), donc la courbe passe par le point de coordonnées \((0;1)\).
  • Ensuite, puisque \(f'(x)=0,5\times f(x)\), alors \(f'(0)=0,5\times f(0)=0,5\), donc l’équation de la tangente à la courbe au point d’abscisse 1 est \(y=0,5 x+1\). Intervient l’approximation de la méthode d’Euler : nous considérons que sur l’intervalle \([0;0,1]\), la courbe de la fonction est confondue avec la tangente, et donc que \(f(0,1)=0,5\times 0,1+1=1,05\). La courbe passe donc par le point de coordonnées \((0,1;1,05)\).
  • Et nous répétons cette étape. Puisque \(f'(x)=0,5\times f(x)\), alors \(f'(0,1)=0,5\times f(0,1)=0,5\times 1,05=0,525\), donc l’équation de la tangente à la courbe au point d’abscisse \(0,1\) est \(y=0,525x+0,9975\). Nous considérons que sur l’intervalle \([0,1;0,2]\), la courbe de la fonction est confondue avec la tangente, donc \(f(0,2)=0,525\times 0,2+0,9975=1,1025\). La courbe passe donc par le point de coordonnées \((0,2;1,1025)\).

Et ainsi de suite.

Les résultats obtenus sur l’intervalle \([0;1]\) sont visibles dans ce fichier (arrondis à \(10^{-4}\)), ou graphiquement sur le tracé suivant (avec en bleu, la solution exacte).

\begin{axis}[% scatter/classes={default={mark=o,draw=black}}, enlargelimits=true, axis x line=center, axis y line=center, xmin=0, ymin=1, ] \addplot[domain=0:1, blue, samples=100]{.5*exp(x)+.5}; \addplot[scatter, scatter src=explicit symbolic] table { x y 0 1 0.1 1.05 0.2 1.1025 0.3 1.1576 0.4 1.2155 0.5 1.2762 0.6 1.3400 0.7 1.4071 0.8 1.4774 0.9 1.5513 1 1.6288 }; \end{axis}

Limites

Comme on peut le voir sur le graphique précédent, l’erreur d’approximation est de plus en plus élevée. Pour améliorer cela, on peut prendre un pas plus petit, mais ce ne sera jamais parfait.

Si je me souviens bien de mes études, un autre problème de la méthode d’Euler est qu’il ne respecte pas la conservation de l’énergie. Si ce n’est pas forcément très grave pour un jeu vidéo, cela peut poser de gros problèmes au moment d’envoyer une fusée dans l’espace. D’autres méthodes sont alors plus adaptées.

Applications

Voici quelques applications mises en œuvre dans ce dépôt.

Chute libre avec rebond

La seconde loi de Newton appliquée à un objet en chute libre donne comme équation \(m\overrightarrow{g}=m\overrightarrow{a}\), soit plus simplement \(\overrightarrow{g}=\overrightarrow{a}\), ainsi, sur chacun des deux axes \(x\) et \(y\) (dans un plan vertical contenant la trajectoire), on a :

\[\begin{split}\left\{\begin{array}{rcl} x''(t) &=& 0\\ y''(t) &=& -9,81\\ \end{array}\right.\end{split}\]

On considère en plus que lorsque l’objet atteint le sol (\(y<0\)), la vitesse est inversée, et amortie.

La mise en œuvre en Python est donc le code suivant.

 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
"""Résolution des équations modélisant une chute libre avec la méthode d'Euler."""

# pylint: disable=invalid-name

import turtle

# Position initiale
px, py = 0, 0
# Vitesse initiale
vx, vy = 5, 40
# Accélération (qui est constance)
ax, ay = 0, -9.81

h = 0.01

while True:
    # Application de la méthose d'Euler
    px, py, vx, vy, ax, ay = (
        px + h * vx,
        py + h * vy,
        vx + h * ax,
        vy + h * ay,
        ax,
        ay,
    )

    # Rebond
    if py < 0:
        py = -py
        vy = -0.9 * vy

    # Tracé
    turtle.goto(px, py)  # pylint: disable=no-member

Ce qui donne l’exemple de trajectoire suivant.

../_images/euler-chute.png

Satellite

La seconde loi de Newton appliquée au satellite d’une étoile donne comme équation \(\overrightarrow{P}=m\overrightarrow{a}\), soit (en ne considérant que l’intensité des forces) \(\frac{MG}{d^2}=a\). Ici, \(M\) et \(G\) (masse de l’étoile, et constante universelle de gravitation) sont constante, mais la distance \(d\) entre les deux objets varie en fonction de la position du satellite. L’équation différentielle est donc (avec \(g\) l’intensité de la gravité divisée par le poids, et \(\alpha\) l’angle formé entre l’axe des abscisses, et la droite étoile-satellite) :

\[\begin{split}\left\{\begin{array}{rcl} x''(t) &=& -g\cos{\alpha}\\ y''(t) &=& -g\sin{\alpha}\\ \end{array}\right.\end{split}\]

La mise en œuvre en Python est donc le code suivant (l’étoile est à l’origine du repère).

 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
"""Résolution des équations modélisant la trajectoire d'un satellite avec la méthode d'Euler."""

# pylint: disable=invalid-name, no-member

import turtle
from math import cos, sin, atan2

# Position initiale (à droite)
px, py = 250, 0
# Vitesse initiale (perpendiculaire à l'axe étoile-satellite)
vx, vy = 0, 40
# Accélération initiale (nulle, mais c'est sans importance)
ax, ay = 0, 0
# Composante M×m×G de la gravité (unité arbitraire)
f = 1000000

h = 0.001

turtle.delay(0.1)
turtle.up()
turtle.goto(px, py)
turtle.down()

while True:
    # Calcul de la gravité
    gravite = f / (px ** 2 + py ** 2)
    # Calcul de l'angle entre l'étoile et son satellite
    angle = atan2(py, px)

    px, py, vx, vy, ax, ay = (
        px + h * vx,
        py + h * vy,
        vx + h * ax,
        vy + h * ay,
        -gravite * cos(angle),
        -gravite * sin(angle),
    )

    turtle.goto(px, py)

Ce qui donne l’exemple de trajectoire suivant.

../_images/euler-satellite.png

Remarquons que si mes souvenirs sont bons, la solution exacte de cette équation différentielle forme une ellipse parfaite : le satillite devrait revenir exactement à son point de départ. Ce n’est pas le cas ici, à cause des approximations de cette méthode.

Ressort amorti

Dans notre système, deux forces sont exercées sur la masse :

  • la force de traction (ou d’extension) du ressort, d’intensité proportionnelle à son allongement (la différence par rapport à la longueur au repos) ;
  • les frottements, d’intensité proportionnelle à la vitesse.

L’équation différentielle à résoudre est donc (où \(k\) et \(f\) sont des constantes, ici arbitraires) :

\[-kx(t)-fx'(t)=mx''(t)\]

Dans ce cas, j’ai choisi de placer le temps sur l’axe des abscisses, et l’allongement sur les ordonnées. La mise en œuvre en Python est donc le code suivant.

 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
"""Résolution des équations modélisant un ressort amorti en utilisant la méthode d'Euler."""

# pylint: disable=invalid-name, no-member

import turtle

# Allongement
x = 300
# Vitesse
v = 0
# Accélération
a = 0
# Temps
t = 0

# Coefficient de frottement
f = 0.5
# Coefficient de la force de rappel
k = 20

h = 0.001

turtle.tracer(8, 1)
size = turtle.screensize()
turtle.up()
turtle.goto(t, x)
turtle.down()

while True:
    t += 0.02

    x, v, a = (x + h * v, v + h * a, -f * v - k * x)

    turtle.goto(t, x)

Ce qui donne l’exemple de trajectoire suivant.

../_images/euler-ressort.png

Système proies-prédateurs

Lotka et Volterra ont proposé, indépendamment, un modèle de l’évolution des populations de proies et prédateurs. Le système d’équations à résoudre est le suivant (voir la signification des constantes dans l’article de Wikipédia cité précédemment).

\[\begin{split}\left\{\begin{array}{rcl} x'(t) &=& x(t) \times \left(\alpha - \beta y(t) \right) \\ y'(t) &=& y(t) \times \left(\delta x(t) - \gamma \right) \end{array}\right.\end{split}\]

La mise en œuvre en Python est donc le code suivant.

 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
48
49
50
51
52
53
"""Résolution des équations de Lotka-Volterra avec la méthode d'Euler."""

# pylint: disable=invalid-name

import turtle

# Définition de trois tortues
t_proie = turtle.Turtle()
t_predateur = turtle.Turtle()
t_phase = turtle.Turtle()

# Populations
n_proie = 50
n_predateur = 20
# Vitesses de variation des populations
v_proie = 0
v_predateur = 0
# Constantes des équations de Lotka Volterra
alpha = 3
beta = 0.1
gamma = 2
delta = 0.1

t = 0
h = 0.001

# Initialisation des tortues
turtle.delay(1)  # pylint: disable=no-member
t_proie.up()
t_proie.goto(t - 200, 5 * n_proie - 200)
t_proie.color("green")
t_proie.down()
t_predateur.up()
t_predateur.goto(t - 200, 5 * n_predateur - 200)
t_predateur.color("red")
t_predateur.down()
t_phase.up()
t_phase.goto(2 * n_proie - 400, 2 * n_predateur)
t_phase.down()

while True:
    t += 0.05

    n_proie, n_predateur, v_proie, v_predateur = (
        n_proie + h * v_proie,
        n_predateur + h * v_predateur,
        n_proie * (alpha - beta * n_predateur),
        n_predateur * (delta * n_proie - gamma),
    )

    t_proie.goto(t - 200, 5 * n_proie - 200)
    t_predateur.goto(t - 200, 5 * n_predateur - 200)
    t_phase.goto(2 * n_proie - 400, 2 * n_predateur)

Ce qui donne l’exemple de trajectoire suivant, où :

  • la courbe verte est celle de la population des proies en fonction du temps ;
  • la courbe rouge est celle de la population des prédateurs en fonction du temps ;
  • la courbe noire est la courbe des points \((x(t);y(t))\).

Nous pouvons observer un cycle.

  • Au départ, les proies sont nombreuses. Les prédateurs ont beaucoup à manger, donc leur population augmente.
  • Puisque la population des prédateurs augmente, ils mangent de plus en plus de proies, dont la population diminue.
  • La population de proies n’est plus assez suffisante pour nourrir les prédateurs, dont la population diminue.
  • Il y a moins de prédateurs, dont la population de proies augmente, ce qui ramène à la population initiale.
../_images/euler-proie-predateur.png