TangenteX.com
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 abordables 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 comment 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 :
Le coefficient directeur de la tangente est ici égal à ( y(x0+
h) - y(x0) ) / ( (x0+ h) - x0 )
ou encore en simplifiant le dénominateur ( y(x0+ h) -
y(x0) ) / h
Vous savez aussi que le coefficient directeur de la tangente en (x0,y0)
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'(x0) = (y(x0+ h) - y(x0))/h
ou encore
y(x0+ h) = y(x0) + h*y'(x0)
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 (x0+ h, y(x0+ h)) en connaissant x0,
y(x0), h et la dérivée de y en x0.
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ériques 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 !
Voyons maintenant comment construire un algorithme qui permette
d'utiliser la méthode d'Euler pour intégrer numériquement une
équation différentielle.
On donne une équation différentielle de la forme y' = f(y) et des
conditions initiales telles que pour y(t = 0) = 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 :
xn = xn+1 + h
yn = yn + y'(xn)*h
et en connaissant x0 et y0.
Commençons par implémenter cet algorithme dans le programme en C EulerCC.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;
}
Tous les codes sources étudiés dans cette page sont téléchargeables dans la bibliotèque de codes de TangenteX.
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
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.
On peut également programmer facilement la méthode d'Euler sur
une TI 89 ou V200. Cherchons un programme 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 emploierons bien sur le même algorithme que 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:
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é !
Voyons maintenant comment employer la méthode d'Euler sur un problème très classique : l'évolution d'une 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 N0 est une constante, la population d'atomes à t0. 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:
Le code source Desintegration.c ci-dessous permet l'intégration de l'équation de désintégration, avec les conditions initiales N0 = 6.0 1020 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;
}
Je vous propose ici le code correspondant en Python dans le script DesintegrationEuler.py, 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)
Voici la courbe obtenue :
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*v2.
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)*v2, où m est la masse du grêlon et k le
coefficient de frottement fluide. Je ne présente 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.
Même principe que pour le problème de désintégration: on adapte le
calcul de la dérivée (ici g-k*v2), les constantes (h,
N) et les conditions initiales.
Il reste à compiler le source grelon.c 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;
}