Methode d'Euler

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

La méthode d'Euler

La méthode d'Euler est une méthode simple 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₀.
Avec cette formule, je peux trouver numériquement une solution à toutes les EDO de premier ordre que l'on rencontre en terminale.
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 exercices de terminale et aussi quelques applications un peu frustres.

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₀.
Il est classique en terminale de traiter cet algorithme à l'aide d'un tableur. Pour ma part, je préfère écrire un petit programme :

//******************************************************************************

//* 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 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é!

Expérience 1 - Désintégration radioactive

Voyons maintenant comment employer la méthode d'Euler sur un problème archi-classique du cours de physique de terminale (programme avant 2012) : l'évolution de la population d'atomes radioactifs.
Vous avez appris dans votre cours que la décroissance de cette population par unité de temps dépendait de deux données : la population à l'instant t, que je note très classiquement N(t) et un paramètre propre au noyau considéré qui s'appelle la constante radioactive (de dimension T^-1 pour information...), que je note λ.
J'écris ces hypothèses sous forme mathématique:

-dN(t)/dt = λ*N(t)

Attention au signe -, qui indique la décroissance de la population. J'obtiens donc l'équation différentielle bien connue:

dN(t)/dt + λ*N(t) = 0

Cette équation possède, comme vous le savez, une solution analytique de la forme N(t) = N₀*e^-^λ^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, ici λ*N(t), soit dx = -λ*x, la variable locale x contenant la valeur de N(t). Je définis évidemment la variable λ 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₀)

Et voilà, je peux compiler et exécuter mon programme. 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 d'atomes 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; // 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 Desintegration.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("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

Expérience 2 - Chute d'une bille avec frottement

Encore un problème archi-classique, toujours avant 2012 : 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 (celui qui en a le courage peut mettre la démonstration sur le forum!). 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. Si c'était le cas, venez en discuter sur le forum.
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; // 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 grelon.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("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