Méthode d'Euler

Il s'agit d'une méthode de résolution numérique des équations différentielles. Par *numérique*, vous devez comprendre qu'on utilise une machine pour trouver une solution approximée.

Une équation différentielle est une équation dans laquelle l'inconnue n'est pas un nombre mais une fonction. On cherche à déterminer cette fonction qui apparaît avec sa (ses) dérivée(s) dans l'équation. C'est un problème important car la plupart des équations de la physique prennent la forme d'une équation différentielle et certaines sont impossibles à résoudre de façon exacte.

Prenons un exemple : $f(t) + 5\cdot f'(t) = 10$ et $f(0)= 0$.

Quelques observations :

• En physique, la variable est souvent le temps $t$.
• En général, au lieu de $f(t)$, on écrit $y$.
• Une dérivée représente un taux d'accroissement : à un instant $t$ donné, la fonction vaut $f(t)$. Si on laisse passer un instant très court $dt$ (notation de physicien) le temps devient $t + dt$ et la fonction $f(t + dt)$. On notera $df = f(t + dt) - f(t)$, la variation de $f$. L'idée de la dérivation est de considérer que, quand $dt$ très petit, $df$ tend à être proportionnel à $dt$. On note $f'(t) = \dfrac{df}{dt}$ ce facteur de proportionnalité, le taux d'accroissement.
• L'équation différentielle est complétée par une condition initiale, ici $f(0) = 0$. En effet, l'équation a une infinité de solutions. Physiquement, chaque solution correspond à une trajectoire possible, une trajectoire permise par les lois de la physique. Mais dans le cadre d'une expérience particulière, une trajectoire a été choisie. La condition initiale permet de sélectionner cette trajectoire.

Dans l'équation donnée en exemple, on ne voit que $f$ et $f'$. C'est un premier ordre. Certaines équations de la physique font apparaître $f$, $f'$ et $f''$. Ce sont des seconds ordre. On peut trouver plus compliqué encore. Dans le cadre de ce TD, on se contentera de premier ordre. Avec les équations du premier ordre, la méthode d'Euler fonctionne très bien. Elle est moins performante avec le second ordre.

L'équation donnée en exemple est simple et on connaît la solution exacte :

$$f(t) = 10\cdot(1 - \exp\left(-\frac{1}{5}t\right)$$

On ne connaît les solutions exactes que de certaines équations différentielles. On aimerait une méthode permettant au moins de connaître une approximation de la solution.

Dans notre exemple, l'équation s'écrit $y + 5 y' = 10$ et $y(0) = 0$.

De plus, $y' = \dfrac{dy}{dt}$. On peut donc réécrire l'équation :

$$y + 5 \frac{dy}{dt} = 10 \Rightarrow dy = \frac{1}{5}(10-y)\cdot dt$$

Pour être mathématiquement exacte, cette équation impose $dt \to 0$. Avec la méthode d'Euler, on prend $dt$ petit mais pas infiniment petit. C'est une approximation. Prenons par exemple $dt = 0,1$.

Nous en arrivons au cœur du fonctionnement de la méthode d'Euler :

• Je connais $y(0)$
• Je sais que dans le temps $dt$, $y$ augmente de la quantité $\frac{1}{5}(10-y)\cdot dt$
• Pour $dt = 0,1$ et $y(0) = 0$ je calcule donc $\frac{1}{5}(10-0)\cdot 0,1 = 0,2$. Cela signifie que $y(0,1) = y(0) + 0,2 = 0,2$.
• Je répète l'opération partant de $t = 0,1$ et $y(0,1) = 0,2$ : $\frac{1}{5}(10-0,2)\cdot 0,1 = 0,196$. Donc $y(0,2) = y(0,1) + 0,196 = 0,396$
• Et on poursuit… jusqu'à atteindre l'instant désiré, par exemple $t = 20$. Pour atteindre $t = 20$ par intervalles de $dt = 0,1$, il nous faudra une boucle répétant 200 fois le calcul.

Nous poursuivons sur notre exemple mais en gardant à l'esprit que l'on voudrait pouvoir résoudre d'autres équations. Nous devons donc mettre notre exemple de sorte qu'on puisse facilement le changer.

On va implémenter l'équation sous sa 2e forme : $dy = \frac{1}{5}(10-y)\cdot dt$.

def exemple(y, dt):
    return 1/5*(10-y)*dt

Cette écriture convient pour l'exemple mais n'est pas assez générale. Voici un autre exemple :

$$4t\cdot y - 8 y' = \sin(2t) \Rightarrow dy = \frac{1}{8}(\sin(2t) - 4t\cdot y)\cdot dt$$

Comme vous le voyez, on a parfois de connaître aussi $t$ pour le calcul. On écrira donc pour l'exemple :

def exemple(y, t, dt):
    return 1/5*(10-y)*dt

Ici la variable t ne sert pas mais on la met tout de même pour rester le plus général possible. Ainsi, si on voulait changer le premier exemple, pour le deuxième exemple on aurait :

# 2e exemple
from math import sin
def exemple(y, t, dt):
    return 1/8*(sin(2t)-4*t*y)*dt

le calcul serait différent mais la signature de la fonction exemple(y, t, dt) serait la même. C'est important si on veut que la suite du programme soit indépendant de l'équation différentielle envisagée. On veut donc une fonction $y, t, dt \mapsto dy$.

On va écrire une fonction euler qui fait la résolution. En entrée de la fonction, on donne :

• la fonction décrivant l'équation différentielle (exemple) qui doit être de forme $y, t, dt \mapsto dy$
• la condition initiale y0 représentant la valeur $y(0)$
• T, le temps total que l'on veut simuler (dans l'exemple, T = 20)
• N, le nombre de pas que l'on veut calculer (dans l'exemple, N = 200)

En sortie la fonction doit renvoyer :

• un tableau contenant tous les instants $t$. Dans l'exemple : [0, 0.1, 0.2, ..., 20]
• un tableau contenant les valeurs correspondantes de $y(t)$. Dans l'exemple [0, 0.2, 0.396, ...]

La fonction suit l'algorithme suivant (en pseudo-code)

DÉBUT
    calcul de dt connaissant T et N
    soit ts = [0] le tableau des temps
    soit ys = [y0] le tableau des y
    RÉPÉTER N FOIS:
        t précédent est dernier temps de ts
        y précédent est dernier y de ys
        t actuel est t précédent + dt
        calcul de dy en utilisant y précédent, t précédent, dt
        y actuel est y précédent + dy
        ajouter t actuel aux ts
        ajouter y actuel aux ys
    FIN RÉPÉTER
    RENVOYER ts et ys
FIN

On peut utiliser matplotlib pour faire le tracer. Puisque euler nous fournit le tableau des temps et des y, c'est très rapide :

import matplotlib.pyplot as plt

def exemple(y, t, dt):
    return 1/5*(10-y)*dt
    
    
def euler(equation, y0, T, N):
    '''
    equation: fonction y, t, dt -> dy
    y0: valeur de y(0)
    T: durée de simulation
    N: nombre d'échantillons
    '''
    # ... votre travail

ts, ys = euler(exemple, 0, 20, 200) # y0 = 0, T = 20, N = 200
plt.plot(ts, ys)
plt.show()

Il est possible d'améliorer l'affichage de la courbe mais c'est un autre travail…