TangenteX.com

Intégrer une équation différentielle par la méthode d'Euler

La méthode d'Euler

Sa mise en oeuvre

La méthode d'Euler en Python

La méthode d'Euler sur une TI 89

Désintégration radioactive

Chute d'une bille avec frottement

Accueil

La méthode d'Euler

La méthode d'Euler est une méthode de résolution d'une équation différentielle ordinaire (EDO) de premier degré. Comme son nom l'indique, elle est due au mathématicien et physicien suisse Euler (1707-1783). Son principe est relativement simple et s'appuie sur des outils tout à fait abordable aux élèves de terminale et même de première, car il s'agit essentiellement de dérivée.
Mais revenons d'abord à la signification d'une équation différentielle de premier ordre. Pour un physicien, une EDO indique l'évolution d'un système physique en fonction de l'accroissement d'un paramètre, très généralement le temps. Voyons sur un graphique comme on peut aborder ça...

Nous avons donc une courbe C de la forme y = f(x). J'ai tracé la tangente à C au point (x0, y0). Vous savez comment lire graphiquement le coefficient directeur d'une droite, ici de la tangente :

• on part d'un point quelconque de la droite d'abscisse x₀,
• on se déplace sur l'axe des x en ajoutant une valeur petite à x₀, valeur que j'appelle ici h,
• et on mesure de combien la position de la nouvelle ordonnée du point x₀+h a variée.

Le coefficient directeur de la tangente est donc ici égal à ( y(x₀+h) - y(x₀) ) / ( (x₀+h) - x₀ ) ou encore en simplifiant le dénominateur ( y(x₀+h) - y(x₀) ) / h
Vous savez aussi que le coefficient directeur de la tangente en (x₀,y₀) est égal à la valeur de la dérivée de y en ce point.
Si h est très petit (en théorie, s'il tend vers 0), je peux écrire l'approximation suivante :

y'(x₀) = (y(x₀+h) -y(x₀))/h

ou encore

y(x₀+h) = y(x₀) + h*y'(x₀)

Attention, il s'agit d'une approximation car h est petit mais ne tend pas vers 0.
J'ai donc une équation qui me permet de calculer la position du point (x₀+h, y(x₀+h)) en connaissant x₀, y(x₀), h et la dérivée de y en x₀.

La méthode d'Euler est simple, mais assez imprécise. On le verra sur un exemple dont on connait la solution analytique. Elle est donc peu employée sous cette forme, sauf pour les applications un peu frustres. Mais n'oublions pas qu'elle est la base de méthodes numérique de résolution des EDO plus élaborées, comme les méthodes Runge Kutta (RK). D'ailleurs, la méthode d'Euler est formellement la méthode RK d'ordre 1 !

Sa mise en oeuvre

Voyons maintenant comment construire un algorithme qui me permette d'utiliser la méthode d'Euler pour intégrer numériquement une équation différentielle.
On me donne une équation différentielle de la forme y' = f(y) et des conditions initiales telles que pour t=0 y(t) = a
Je vais choisir un pas d'intégration, que j'appelle h, aussi petit que possible. Et là, il va falloir faire un compromis entre précision et temps de calcul. Plus h est petit et plus ma solution sera précise. Mais plus le temps de calcul sera long! Et de plus, la méthode souffre d'autres défauts qui font que de toute manière le résultat ne sera pas très précis. On verra ce problème dans les exemples.
Je vais donc construire une courbe solution de mon équation différentielle, pas à pas, en remarquant que :

x_n = x_n+1 + h

y_n = y_n + y'(x_n)*h

et en connaissant x₀ et y₀.

Commençons par implémenter cet algorithme dans un petit programme en C :

//******************************************************************************
//* Programme de démonstration d'utilisation de la méthode d'Euler
//* pour résoudre une EDO de premier ordre.
//*
//* Dominique Lefebvre - TangenteX.com octobre 2008
//******************************************************************************
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//******************************************************************************
//* Calcul de la valeur de la dérivée en x
//* l'équation différentielle est y' = 1/x²
//******************************************************************************
double Derivee(double x)
  double dx;
  dx = 1/(x*x);
return dx;
}
//******************************************************************************
//* Programme principal
//******************************************************************************
int main(int argc, char *argv[])
{
  //* Declaration des variables
  int i, N;
  double h;
  double *x, *y;
  FILE *fp;
  //* Déclaration des constantes de calcul
  h = 0.01; // pas pour le calcul
  N = 500;
  //* Allocation de la mémoire pour les tableaux de calcul
  x = (double *)malloc((unsigned)(N+1)*sizeof(double));
  y = (double *)malloc((unsigned)(N+1)*sizeof(double));
  //* Determination des conditions initiales x0= 1 y0 = 0
  x[0] = 1.0; // abscisse initiale
  y[0] = 0.0; // ordonnée initiale
  //* Calcul
  for (i=0; i<N; i++)
   {
     y[i+1] = y[i] + Derivee(x[i])*h; // on applique le schéma d'Euler
     x[i+1] = x[i] + h;
   }
  // Sauvegarde des résultats dans un fichier texte pour tracer par GnuPlot
  // Ce fichier se nomme EulerC.dat. Il est cree dans le repertoire courant
  // (celui dans lequel est lance le programme)
  // S'il n'existe pas, il est cree. S'il existe, son contenu est ecrase.
  fp=fopen("EulerC.dat","w+");
  for(i=0;i<N; i++)
    fprintf(fp,"%f %f\n",x[i], y[i]);
  fclose(fp);
  //* Liberation de la memoire et sortie
  free(x);
  free(y);
  system("PAUSE");
  return 0;
}

Le code source est téléchargeable : EulerCC.c. Vous pouvez le compiler et l'exécuter avec DevC++ ou n'importe quel environnement de développement C (Visual studio, Code::Blocks ou autre).

La méthode d'Euler en Python

La méthode d'Euler peut être également implémentée très simplement avec Python. Commençons par définir la fonction Euler(), qui implémente l'algorithme :

def Euler(fonction,y0,t):
    """ Intégration d'une EDO du premier ordre par la méthode d'Euler
        fonction : EDO à intégrer
        y0 : valeur initiale
        t : vecteur variable d'intégration
    """
    # détermination des paramètres d'intégration nombre de pas (n) et
    # intervalle entre deux pas (h)
    n = len(t)
    h = t[1]-t[0]
    # initialisation du vecteur de retour
    y = zeros(n)
    y[0] = y0  
    # calcul du vecteur de retour
    for i in range(n-1):
        y[i+1] = y[i] + fonction(y[i], t[i]) * h
    return y

Cette fonction est disponible dans le package Python TxEDO.py. Ce package contient aussi une fonction d'intégration des EDO par RK4.

Puis reprenons l'exemple ci-dessus, en le codant cette fois en Python plutôt qu'en C. Voici le code correspondant :

# importation des librairies
from TxEDO import Euler, VecteurTemps
from matplotlib.pyplot import figure,plot,xlim,xlabel,ylabel,grid,title

# définition de l'EDO à intégrer
def EDO(x,t):
    return (1.0/x**2)
   
# définition du vecteur temps de l'expérience
t0 = 0.0
tmax = 500.0
dt = 0.01
t = VecteurTemps(t0, tmax, dt)

# définition de la condition initiale
x0 = 1.0

# solution par la méthode d'Euler
y = Euler(EDO,x0,t)

# Affichage de la solution
figure(1)
grid(True)
plot(t, y,'blue')
xlim(t0,tmax)
xlabel('$temps$', fontsize = 15)
ylabel('$y$', fontsize = 15)
title('Integration 1/x^2', fontsize = 15)

Vous constatez que ce code est plus "léger" que le code C ci-dessus, d'autant qu'il intégre l'affichage de la courbe résultante, mais il respecte la même structure. Il est disponible dans le script ExempleEuler.py.

La méthode d'Euler sur une TI 89

On peut également programmer facilement la méthode d'Euler sur une TI 89 ou V200. Cherchons un programme simple qui permette de tracer une solution de l'équation différentielle y'(t) = f(y,t) avec une condition initiale y(t0) = Cts, ce que l'on appelle un problème de Cauchy..
Nous employerons bien sur le même algorithme que décrit ci-dessus, mais le langage de programmation de la TI 89 (ou V200) rendra l'écriture du programme un peu plus simple! Le programme se réduit à la liste d'instructions suivantes:

:euler()
:Prgm
:FnOff
:Local t,y,p
:Input "t0: ",t
:Input "y0: ",y
:Input "pas calcul:",p
:Lbl BP1
:For i,0,20
:t+p -> u
:y + y1(t)*p -> v
:Line t,y,u,v
:u -> t
:v -> y
:EndFor
:Pause
:Goto BP1
:EndPrgm


Ce programme appelle peu de commentaires:

• les caractères -> désignent la commande STO de la TI 89
• il vous demande les conditions initiales, soit la valeur y0 de y pour t = t0 et le pas de calcul
• le programme affiche 20 pas de calcul puis attend que vous appuyez sur ENTER pour afficher les 20 suivants.
• vous devez entrer l'équation y' = f(y,t) sur la ligne y1 de l'éditeur d'équation
• vous devez également définir les coordonnées d'affichage dans WINDOW

Je suppose que vous savez entrer un programme dans votre calculatrice, soit directement soit par un éditeur sur votre PC. Sinon contactez-moi.
Comme vous le constatez, l'implémentation de la méthode d'Euler sur une TI n'est pas très compliqué!

Désintégration radioactive

Voyons maintenant comment employer la méthode d'Euler sur un problème très classique : l'évolution de la population de noyaux radioactifs.

Je ne reviens pas sur la définition de la désintégration radioactive dont vous pourrez trouver une description sur la page de TangenteX qui lui est consacrée. Je retiendrai simplement que l'équation différentielle modélisant ce phénomène est :

\( \dfrac{dN}{dt} + \lambda N = 0 \)

Cette équation possède, comme vous le savez, une solution analytique de la forme \( N(t) = N_0 e^{-\lambda t} \), dans laquelle N₀ est une constante, la population d'atomes à t₀ en l'occurrence. Tant mieux, cela vous permettra de vérifier les résultats obtenus en utilisant la méthode d'Euler pour résoudre numériquement cette équation différentielle.
Je vais donc modifier le programme exemple donné ci-dessus pour l'adapter à ce cas particulier. Et cette adaptation est très simple:

• dans la routine Derivee, j'indique la formule de la dérivée de mon problème, la variable locale x contenant la valeur de N(t). Je définis évidemment la variable \( \lambda \) en fonction de sa valeur.
• dans le corps du programme principal, j'ajuste les valeurs de h (pas de calcul) et de N (durée de l'expérience). La valeur de h doit rester petite par rapport à la durée de l'expérience. Un pas de 0,1 s dans ce cas est correct. Mais vous pouvez le faire varier pour examiner l'influence de la valeur de h sur la précision du résultat. La valeur de N dépend de vos besoins.
• toujours dans le programme principal, je définis les conditions initiales. Elles sont au nombre de deux dans notre problème:
  • le temps, porté par la variable x. On choisit t₀ = 0.
  • la population initiale, portée par la variable y (à t₀, N(t) = y = N₀)

Son implémentation en C

Le code C ci-dessous permet l'intégration de l'équation de désintégration, avec les conditions initiales N₀ = 6.0 10²⁰ noyaux et la constante de désintégration égale à 1.8 10^-3. Les résultats sont tracés dans un fichier (ici desintegration.dat) que vous visualiserez en utilisant GnuPlot.

//******************************************************************************
//* Programme de calcul de l'évolution d'une population de noyaux radioactifs
//* en utilisant la méthode d'Euler.
//*
//* Dominique Lefebvre - TangenteX.com octobre 2008
//******************************************************************************
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//******************************************************************************
//* Calcul de la valeur de la dérivée en x
//* l'équation différentielle est y' = -lambda*x
//******************************************************************************
double Derivee(double x)
{
double dx;
double lambda = 1.8e-3;
dx = -lambda*x;
return dx;
}
//******************************************************************************
//* Programme principal
//******************************************************************************
int main(int argc, char *argv[])
{
  // Declaration des variables
  int i, N;
  double h;
  double *x, *y;
  FILE *fp;
  //* Déclaration des constantes de calcul
  h = 0.1; // pas pour le calcul
  N = 100;
  // Allocation de la mémoire pour les tableaux de calcul
  x = (double *)malloc((unsigned)(N+1)*sizeof(double));
  y = (double *)malloc((unsigned)(N+1)*sizeof(double));
  // Determination des conditions initiales x0= 1 y0 = 0
  x[0] = 0.0; // t0 = 0
  y[0] = 6.00e20; // N0 :nombre d'atomes à t0
  //* Calcul
  for (i=0; i<N; i++)
   {
     y[i+1] = y[i] + Derivee(y[i])*h;
     x[i+1] = x[i] + h;
   }
  // Sauvegarde des résultats
  fp=fopen("Desintegration.dat","w+");
  for(i=0;i<N;i++)
    fprintf(fp,"%f %f\n",x[i], y[i]);
  fclose(fp);
  //* Liberation de la memoire et sortie
  free(x);
  free(y);
  system("PAUSE");
  return 0;
}

Le code source de ce programme est disponible dans le fichier Desintegration.c

Son implémentation en Python

Je vous propose ici le code correspondant en Python, avec les mêmes paramètres et conditions initiales :

# importation des libraires
from TxEDO import Euler, VecteurTemps
from matplotlib.pyplot import figure,plot,xlim,xlabel,ylabel,grid,title
   
# définition de l'EDO à integrer   
def Desintegration(x,t):
    x_point = -Lambda*x
    return x_point
   
# définition du vecteur temps de l'expérience
t0 = 0.0
tmax = 100.0
dt = 0.1
t = VecteurTemps(t0, tmax, dt)

# définition des paramètres de l'expérience
Lambda = 1.8e-3     # constante de désintégration
N0 = 6.00e20        # nombre de noyaux initial

# solution par la méthode d'Euler
N_Euler = Euler(Desintegration,N0,t)

# Affichage des deux solutions
figure(1)
grid(True)
plot(t, N_Euler,'blue')
xlim(t0,tmax)
xlabel('$temps$', fontsize = 15)
ylabel('$Population$', fontsize = 15)
title('Desintegration', fontsize = 15)

Le script est disponible dans DesintegrationEuler.py. Voici la courbe obtenue :

Chute d'une bille avec frottement

Encore un problème très classique : la chute d'une bille (un grêlon par exemple) avec une force de frottement fluide, pour laquelle on cherche la vitesse limite.
Dans le cas d'un grêlon, on montrera facilement que la poussée d'Archimède est complètement négligeable. D'autre part, la viscosité de l'air et les vitesses en jeu poussent à envisager une force de frottement fluide en k*v².
En appliquant la seconde loi de Newton dans le référentiel adéquat, on obtient donc l'équation différentielle : dv/dt = g -(k/m)*v², où m est la masse du grêlon et k le coefficient de frottement fluide. Je ne présenta pas g.. Là encore, j'imagine que l'établissement de cette équation différentielle ne pose pas de problème.
Le but du jeu est de trouver la vitesse limite numériquement (sa détermination analytique vous permettra de contrôler vos résultats numériques). On va donc calculer la courbe de variation de la vitesse en fonction du temps. L'asymptote horizontale donnera la valeur de la vitesse limite.
Pour cela, on construit le programme à partir du programme exemple donné ci-dessus. Même principe que pour le problème de désintégration: on adapte le calcul de la dérivée (ici g-k*v²), les constantes (h, N) et les conditions initiales.
Il reste à compiler et exécuter le programme. Les résultats sont tracés dans un fichier (ici grelon.dat) que vous visualiserez en utilisant GnuPlot.

//******************************************************************************
//* Programme de calcul de Calcul de la chute d'un grélon
//* en utilisant la méthode d'Euler.
//*
//* Dominique Lefebvre - TangenteX.com octobre 2008
//******************************************************************************
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//******************************************************************************
//* Calcul de la valeur de la dérivée en x
//* l'équation différentielle est y' = k*x²
//******************************************************************************
double Derivee(double x)
  //* definition des constantes
  g = 9.81;
  k = 1.56e-2;
  dx = g-k*(x*x);
return dx;
}
//******************************************************************************
//* Programme principal
//******************************************************************************
int main(int argc, char *argv[])
{
  // Declaration des variables
  int i, N;
  double h;
  double *x, *y;
  FILE *fp;
  //* Déclaration des constantes de calcul
  h = 0.1; // pas pour le calcul
  N = 1000;
  // Allocation de la mémoire pour les tableaux de calcul
  x = (double *)malloc((unsigned)(N+1)*sizeof(double));
  y = (double *)malloc((unsigned)(N+1)*sizeof(double));
  // Determination des conditions initiales x0= 1 y0 = 0
  x[0] = 0.0; // abscisse initiale
  y[0] = 0.0; // ordonnée initiale
  //* Calcul
  for (i=0; i<N; i++)
    {
      y[i+1] = y[i] + Derivee(y[i])*h;
      x[i+1] = x[i] + h;
    }
  // Sauvegarde des résultats
  fp=fopen("grelon.dat","w+");
  for(i=0;i<N;i++)
    fprintf(fp,"%f %f\n",x[i], y[i]);
  fclose(fp);
  //* Liberation de la memoire et sortie
  free(x);
  free(y);
  system("PAUSE");
  return 0;
}

Le code source de ce programme est disponible dans le fichier Grelon.c

Contenu et design par Dominique Lefebvre - www.tangenteX.com janvier 2019 Contact :