Intégration des équations différentielles : méthode d'Euler

3.a. Définition

L'intégration numérique de l'équation $\mathrm{(6)}$ consiste à calculer des valeurs approchées de y sur un intervalle [0,T]. Pour cela, on divise cet intervalle en N sous-intervalles égaux de longueur h=T/N et on définit les instants :

$$t_n=nh\text{(9)}$$

où l'entier n varie de 0 à N. h est le pas de temps. La première valeur Y₀ est la condition initiale.

Dans la méthode d'Euler explicite, la valeur approchée à l'instant t_n+1 est obtenue à partir de la précédente par

$$Y_{n+1}=Y_n+hf(Y_n,t_n)\text{(10)}$$

La méthode d'Euler repose sur l'observation suivante. Supposons que la valeur exacte de Y à l'instant t_n soit connue. La valeur de Y à l'instant t_n+1=t_n+h est donnée par le développement de Taylor :

$$Y(t_n+h)=Y\left(t_n\right)+Y\text{'}\left(t_n\right)h+Y\text{'}\text{'}\left(t_n\right)\frac{1}{2}{h}^2+\dots \text{(11)}$$

La méthode d'Euler consiste à négliger dans ce développement tous les termes à partir du terme d'ordre 2 en h :

$$Y\left(t_{n+1}\right)\simeq Y\left(t_n\right)+hY\text{'}\left(t_n\right)\text{(12)}$$

L'erreur commise par cette approximation (la somme des termes à partir du rang 2) est appelée erreur de troncature (ET), ou erreur locale. Pour la méthode d'Euler, l'ET est O(h²) car si h est assez petit il existe une constante K telle que $\left|ET\right|\le K{h}^2$.

Cela conduit à la relation de récurrence $\mathrm{(10)}$. Il faut remarquer que Y_n est déjà une valeur approchée à l'instant t_n et non pas la valeur exacte. La dérivée est évaluée avec cette valeur approchée; elle est donc en général légèrement différente de la dérivée exacte à l'instant t_n.

Le système différentiel est défini dans une fonction python de la forme systeme(Y,t) où Y est un tableau contenant les variables (x et v dans l'exemple). Par convention, on utilise des lettres majuscules pour les tableaux. La fonction renvoie un tableau numpy.ndarray contenant les dérivées des variables par rapport au temps.

Voici par exemple la fonction définissant le système pour l'oscillateur harmonique (on pose ω=2π) :

import numpy

w2 = (2*numpy.pi)**2
def oscillateur(Y,t):
    return numpy.array([Y[1],-w2*Y[0]])
		      

Par convention, nous avons placé en premier paramètre le tableau Y, en second le temps t. Cette convention est aussi utilisée par la fonction scipy.integrate.odeint. En revanche, la fonction scipy.integrate.solve_ivp utilise la convention inverse (temps en premier, tableau des variables en second).

La fonction pas_euler(systeme,h,tn,Yn) effectue un pas élémentaire de la méthode d'Euler, c'est-à-dire le calcul de Y_n+1 à partir de Y_n. Ses arguments sont le système différentiel, le pas de temps h, le temps tn et le tableau Yn. La fonction renvoie Y_n+1.

def pas_euler(systeme,h,tn,Yn):
    deriv = systeme(Yn,tn)
    return Yn+h*deriv
		       

La fonction euler(systeme,Yi,T,h) effectue l'intégration numérique sur l'intervalle [0,T]. La condition initiale est Yi. Le pas de temps est h. La fonction stocke les valeurs de Y sous forme d'une liste (construite comme une pile). Une liste contenant le temps est aussi générée. Ces listes sont converties en tableau numpy avant d'être renvoyées. Voici par exemple la structure du tableau Y pour un système à 3 variables :

$$\begin{array}{ccc}\hfill Y_0\left[0\right]\hfill & \hfill Y_0\left[1\right]\hfill & \hfill Y_0\left[2\right]\hfill \\ \hfill Y_1\left[0\right]\hfill & \hfill Y_1\left[1\right]\hfill & \hfill Y_1\left[2\right]\hfill \\ \hfill Y_2\left[0\right]\hfill & \hfill Y_2\left[1\right]\hfill & \hfill Y_2\left[2\right]\hfill \end{array}\text{(13)}$$

L'indice ajouté entre crochets permet de repérer la variable. Ainsi Y_n[0] désigne la valeur de la variable d'indice 0 à l'instant n.

def euler(systeme,Yi,T,h):
    Y = Yi
    t = 0.0
    liste_y = [Y]
    liste_t = [t]
    while t<T:
        Y = pas_euler(systeme,h,t,Y)
        t += h
        liste_y.append(Y)
        liste_t.append(t)
    return (numpy.array(liste_t),numpy.array(liste_y))
		       

Cette fonction utilise deux listes remplies au fur et à mesure avec append. Dans le cas présent, la longueur du tableau final est connue a priori (égale à T/h), mais ce n'est pas toujours le cas, par exemple si on met on place une adaptation du pas de temps, ou si le calcul doit s'arrêter lorsqu'une condition sur une des variables est remplie. Le résultat est renvoyé sous la forme de deux tableaux numpy.ndarray, ce qui facilite les manipulations ultérieures, par exemple pour le tracé des courbes. Il est possible d'écrire une version qui travaille entièrement sur ce type de tableaux :

def euler_bis(systeme,Yi,T,h):
    Y = Yi
    t = 0.0
    tableau_y = numpy.array([Y])
    tableau_t = numpy.array([t])
    while t<T:
        Y = pas_euler(systeme,h,t,Y)
        t += h
        tableau_y = numpy.append(tableau_y,[Y],axis=0)
        tableau_t = numpy.append(tableau_t,t)
    return (tableau_t,tableau_y)
		       

Voici le résultat avec l'oscillateur harmonique :

from matplotlib.pyplot import *
T = 5.0
h = 1.0e-3
Yi = [1.0,0]
(t,tab_y) = euler(oscillateur,Yi,T,h)
x = tab_y[:,0]
figure(figsize=(8,4))
plot(t,x)
xlabel('t')
ylabel('y[0]')
axis([0,T,-2,2])
title('h=0.001')
grid()
		       

figA.pdf

3.b. Erreurs locale et globale

L'approximation $\mathrm{(12)}$ consiste à négliger les termes d'ordre supérieur ou égal à 2 dans le développement de Taylor. L'erreur commise par cette approximation est appelée erreur locale de troncature. Lorsque h est assez petit, cette erreur est en valeur absolue inférieure à une constante multipliée par h². Si le pas de temps h est divisé par deux, il y a une réduction d'un facteur 4 de l'erreur de troncature.

L'erreur globale (ou simplement l'erreur) à l'instant t_n est l'écart entre la solution approchée et la solution exacte à cet instant :

$$e_n=y\left(t_n\right)-y_n\text{(14)}$$

L'erreur globale résulte de l'accumulation des erreurs locales commises sur les n pas précédents (ce qui ne signifie pas qu'elle soit toujours croissante). La figure suivante montre un exemple avec la solution exacte en trait plein et les valeurs approchées aux instants t_n.

Figure pleine page

Dans le cas de l'oscillateur harmonique, l'erreur peut être calculée puisque la solution exacte est connue. On trace donc l'erreur en fonction du temps :

erreur = x-numpy.cos(2*numpy.pi*t)
figure(figsize=(8,4))
plot(t,erreur,label='h=0.001')
xlabel('t')
ylabel('e')
grid()
             

figB.pdf

L'erreur présente ici une variation non monotone. Elle s'annule même périodiquement. Globalement, on constate néanmoins une augmentation de l'erreur.

3.c. Convergence

Une méthode numérique doit être convergente, c'est-à-dire que l'erreur globale doit, pour tout n, tendre vers zéro lorsque le pas h tend vers zéro. La méthode d'Euler est effectivement convergente.

La convergence nous garantit des valeurs approchées de plus en plus proches des valeurs exactes lorsqu'on réduit le pas de temps h.

Une question importante est celle de la vitesse de convergence, c'est-à-dire comment l'erreur globale évolue avec le pas h.

Divisons le pas par deux pour voir l'évolution de l'erreur :

h = h/2
(t2,tab_y) = euler(oscillateur,Yi,T,h)
yt = tab_y.transpose()
x = yt[0]
erreur2 = x-numpy.cos(2*numpy.pi*t2)
plot(t2,erreur2,label='h=0.0005')
h = h/2
(t3,tab_y) = euler(oscillateur,Yi,T,h)
x = tab_y[:,0]
erreur3 = x.copy()
for i in range(x.size):
    erreur3[i] = x[i]-math.cos(2*math.pi*t3[i])
plot(t3,erreur3,label='h=0.00025')
legend(loc='upper right')
		          

figC.pdf

L'erreur semble divisée par deux lorsque le pas est divisé par deux. Pour le confirmer, nous allons calculer l'erreur pour une valeur particulière de t, là ou elle présente un maximum relatif. Pour cela, il est inutile de stocker tous les points antérieurs. On écrit une fonction erreur_euler(systeme,solution,yi,T,h) qui calcule la valeur absolue de l'erreur à l'instant final T, pour la première variable. La solution exacte solution doit être fournie.

def erreur_euler(systeme,solution,Yi,T,h):
        Y = Yi
        t = 0.0
        while t<T:
            Y = pas_euler(systeme,h,t,Y)
            t += h
        return abs(Y[0]-solution(t))
		          

On trace l'erreur en t=5 en fonction de h, en échelle logarithmique. Pour cela, on génère une liste de valeurs de h et une liste contenant les erreurs correspondantes.

liste_h = []
liste_e = []
T = 5.0
h = 0.001
def solution(t):
    return math.cos(2*math.pi*t)
for i in range(10):
    liste_h.append(h)
    liste_e.append(erreur_euler(oscillateur,solution,Yi,T,h))
    h = h/2
figure(figsize=(8,6))
plot(liste_h,liste_e,'o')
xlabel('h')
ylabel('|e|')
xscale('log')
yscale('log')
grid()
		          

figD.pdf

On obtient en échelle logarithmique une pente de 1, ce qui confirme que l'erreur globale est proportionnelle à la puissance 1 du pas h. La méthode d'Euler est donc d'ordre 1. Pour réduire l'erreur d'un facteur 10, il faut aussi réduire le pas d'un facteur 10, c'est-à-dire faire 10 fois plus de calculs. Une méthode d'ordre 1 est donc à éviter si l'on souhaite faire des calculs très précis. On utilisera d'autres méthodes d'ordre plus élevé, par exemple la méthode du point milieu (d'ordre 2), une méthode de Runge-Kutta d'ordre 4, ou la méthode de Bulirsch-Stoer.

3.d. Stabilité

On s'intéresse à présent à l'évolution de l'erreur avec le temps T, à h constant.

T = 100.0
h = 1.0e-3
(t,tab_y) = euler(oscillateur,Yi,T,h)
x = tab_y[:,0]
erreur = numpy.absolute(x-numpy.cos(2*numpy.pi*t))
figure(figsize=(8,5))
plot(t,erreur)
xlabel('t')
ylabel('|e|')
yscale('log')
xscale('log')
grid()
		         

figE.pdf

Bien que l'erreur soit très faible au début, elle croît très rapidement. La méthode d'Euler est dite instable pour l'oscillateur harmonique car l'erreur est non bornée.

Voyons le tracé dans l'espace des phases :

v = tab_y[:,1]
figure(figsize=(6,6))
plot(x,v)
xlabel('x')
ylabel('v')
			

figF.pdf

On obtient une oscillation dont l'amplitude augmente sans fin, au lieu du mouvement périodique de la solution exacte. L'instabilité conduit dans ce cas à une solution physiquement inacceptable (l'énergie n'est pas conservée). On peut effectivement démontrer que la méthode d'Euler est instable pour l'oscillateur harmonique (voir paragraphe 6). L'expérience montre qu'elle l'est aussi pour des systèmes non linéaires en général. Pour la simulation des systèmes physiques sur des temps longs, c'est un inconvénient bien plus grave que la faible précision. Lorsqu'on simule un système physique complexe, on recherche plus la vraisemblance de la solution que sa précision. Dans le cas des systèmes chaotiques, la recherche de la précision est d'ailleurs vaine, en raison de l'extrême sensibilité aux conditions initiales.

Des méthodes plus précises comme les méthodes de Runge-Kutta (explicites) sont également instables pour les systèmes linéaires. La précision et la stabilité sont deux propriétés indépendantes.

Nous allons voir comment une modification de la méthode d'Euler permet de la rendre stable.

4.a. Définition

Cette variante de la méthode d'Euler s'applique aux équations différentielles des systèmes mécaniques conservatifs (systèmes hamiltoniens). Elle consiste à traiter séparément les positions et les vitesses. Voyons comment cette méthode est définie dans le cas particulier d'une variable x et de sa vitesse v (la généralisation à un système quelconque est immédiate).

La première étape est le calcul de x_n+1 comme pour la méthode d'Euler :

$$x_{n+1}=x_n+h{v}_n\text{(15)}$$

Avant de calculer la vitesse à l'instant t_n+1, on calcule sa dérivée (c.a.d. l'accélération) en utilisant la valeur x_n+1 de la position. Si on note a(x) la fonction qui permet de calculer l'accélération en fonction de la position, on calcule donc :

$$a_{n+1}=a\left(x_{n+1}\right)\text{(16)}$$

Cela est faisable pour les systèmes conservatifs car la force ne dépend pas de la vitesse. Enfin on calcule la vitesse à l'instant t_n+1 en utilisant cette accélération :

$$v_{n+1}=v_n+h{a}_{n+1}\text{(17)}$$

Dans la méthode d'Euler classique, l'accélération utilisée pour calculer v_n+1 est évaluée avec x_n et non pas x_n+1.

Cette méthode est asymétrique car elle traite de manière asymétrique les variables positions et les variables vitesses. Elle appartient à la classe des méthodes symplectiques (une propriété qui permet de conserver l'énergie). On la désignera sous le nom de méthode d'Euler-A (car il existe aussi une méthode Euler-B, où l'asymétrie est inversée).

Nous allons voir que ce changement modifie radicalement la stabilité de la méthode.

Les équations sont définies par une fonction accel(X). La liste X contient les positions. Voici la fonction pour l'oscillateur harmonique :

def accel(X):
    return numpy.array([-w2*X[0]])
	                 

La fonction pas_eulerA(accel,h,t,Xn,Vn) effectue un pas élémentaire de la méthode Euler-A. La fonction renvoie X_n+1 et V_n+1.

def pas_eulerA(accel,h,t,Xn,Vn):
    X = Xn+h*Vn
    A = accel(X)
    V = Vn+h*A
    return (X,V)
	                   

La fonction eulerA(accel,Xi,Vi,T,h) effectue l'intégration numérique sur l'intervalle [0,T] avec la condition initiale définie par les positions Xi et Vi et un pas de temps h.

def eulerA(accel,Xi,Vi,T,h):
    X = numpy.array(Xi)
    V = numpy.array(Vi)
    t = 0.0
    liste_t = [t]
    liste_x = [X]
    liste_v = [V]
    while t<T:
        (X,V) = pas_eulerA(accel,h,t,X,V)
        t += h
        liste_t.append(t)
        liste_x.append(X)
        liste_v.append(V)
    return (numpy.array(liste_t),numpy.array(liste_x),numpy.array(liste_v))
	                     

4.b. Stabilité

Voyons l'application à l'oscillateur harmonique, en reprenant la valeur de h qui donnait une divergence très visible avec la méthode d'Euler.

T = 1000.0
h = 1.0e-3
Xi = [1]
Vi = [0]
(t,tab_x,tab_v) = eulerA(accel,Xi,Vi,T,h)
figure(figsize=(6,6))
plot(tab_x[:,0],tab_v[:,0])
xlabel('x')
ylabel('v')
title('h=0.001')

	                      

figG.pdf

On obtient bien un mouvement périodique, même avec une durée T très longue.

Voyons le tracé de l'erreur :

x = tab_x[:,0]
erreur = x-numpy.cos(2*numpy.pi*t)
figure(figsize=(8,4))
plot(t,erreur)
xlabel('t')
ylabel('e')
title('h=0.001')
grid()
	                      

figH.pdf

L'erreur garde des valeurs très faibles même sur des temps très longs, là où la méthode d'Euler donnait une divergence catastrophique de l'erreur. On peut démontrer que cette méthode est stable pour l'oscillateur harmonique si h est inférieur à 2/ω (voir paragraphe 6 ci-dessous).

7.a. Principe

L'étude de la stabilité d'un schéma numérique peut se faire de manière générale dans le cas d'un système différentiel linéaire ([3]), qui s'écrit sous la forme :

$$\frac{dY}{dt}=MY\text{(21)}$$

On suppose que la matrice M est diagonalisable. On peut alors effectuer un changement de base qui conduit à résoudre le système suivant :

$$\frac{d\tilde{Y}}{dt}=\tilde{M}\tilde{Y}\text{(22)}$$

où la matrice $\tilde{M}$ est diagonale. On obtient ainsi des équations différentielles découplées, de la forme :

$$\frac{d\tilde{y}}{dt}=q\tilde{y}\text{(23)}$$

où q est une valeur propre de M. La solution de cette équation est :

$$\tilde{y}\left(t\right)=\tilde{y}\left(0\right)exp\left(qt\right)\text{(24)}$$

Les fonctions y(t) solutions du système différentiel sont des combinaisons linéaires de ces fonctions.

Pour l'étude de la stabilité, on s'intéresse aux solutions bornées. Il faut pour cela que les valeurs propres q aient toutes une partie réelle négative ou nulle. Par exemple, pour l'oscillateur harmonique les valeurs propres sont q₁=iω et q₂=-iω et la solution est bornée.

Pour un système dont la solution est bornée, le schéma numérique est stable s'il donne également une solution bornée.

Considérons alors le schéma numérique sous forme matricielle :

$$Y_{n+1}=P{Y}_n\text{(25)}$$

La matrice P est le propagateur. Par exemple pour la méthode d'Euler appliquée à l'oscillateur harmonique :

$$P=\left(\begin{array}{cc}\hfill 1\hfill & \hfill h\hfill \\ \hfill -{\omega }^{2}h\hfill & \hfill 1\hfill \end{array}\right)\text{(26)}$$

La diagonalisation du propagateur (lorsqu'elle est possible), conduit par un changement de base à :

$${\tilde{Y}}_{n+1}=\tilde{P}{\tilde{Y}}_n\text{(27)}$$

où la matrice $\tilde{P}$ est diagonale. On obtient ainsi des équations découplées de la forme :

$${\tilde{y}}_{n+1}=p{\tilde{y}}_n\text{(28)}$$

où p est une valeur propre du propagateur. La solution du schéma numérique est donc :

$${\tilde{y}}_n=p^n{\tilde{y}}_0\text{(29)}$$

La valeur propre p est en général un nombre complexe. Pour que la solution ci-dessus soit bornée, il faut et il suffit que son module soit inférieur ou égal à 1.

On arrive ainsi à la règle suivante : un schéma numérique linéaire est stable pour un système différentiel linéaire dont la solution est bornée si et seulement si toutes les valeurs propres du propagateur ont un module inférieur ou égal à 1.

7.b. Méthode d'Euler et oscillateur harmonique

Appliquons la règle précédente à la méthode d'Euler et à l'oscillateur harmonique. Le propagateur est :

$$P=\left(\begin{array}{cc}\hfill 1\hfill & \hfill h\hfill \\ \hfill -{\omega }^{2}h\hfill & \hfill 1\hfill \end{array}\right)\text{(30)}$$

Ses valeurs propres sont :

$$\begin{array}{ccc}\hfill & p_1=1+ih\omega \hfill & \hfill \text{(31)}\\ \hfill & p_2=1-ih\omega \hfill & \hfill \text{(32)}\end{array}$$

Elles ont un module toujours supérieur à 1. La méthode d'Euler conduit donc à une suite non bornée pour l'oscillateur harmonique, alors que la solution du système différentiel est bornée (fonction cosinus). On en déduit que la méthode d'Euler est instable pour l'oscillateur harmonique.

7.c. Méthode d'Euler asymétrique et oscillateur harmonique

Écrivons explicitement le schéma de la méthode d'Euler asymétrique (Euler-A) pour l'oscillateur harmonique :

$$\begin{array}{ccc}\hfill & x_{n+1}=x_n+h{v}_n\hfill & \hfill \text{(33)}\\ \hfill & v_{n+1}=v_n+h(-{\omega }^2{x}_{n+1})=(1-h^2{\omega }^2)v_n-h{\omega }^2{x}_n\hfill & \hfill \text{(34)}\end{array}$$

Le propagateur est :

$$P=\left(\begin{array}{cc}\hfill 1\hfill & \hfill h\hfill \\ \hfill -{\omega }^{2}h\hfill & \hfill 1-{\omega }^2{h}^2\hfill \end{array}\right)\text{(35)}$$

Les valeurs propres sont solutions de l'équation :

$$p^2+(h^2{\omega }^2-2)p+1=0\text{(36)}$$

dont le discriminant est :

$$\Delta =h^2{\omega }^2(h^2{\omega }^2-4)\text{(37)}$$

Si hω<2 on obient deux racines complexes conjuguées dont le module est égal à 1. Dans ce cas, le schéma est stable. Si hω>2, on obtient deux racines réelles dont le module au carré est :

$$|p{|}^2=1+h^2{\omega }^2(\frac{h^2{\omega }^2}{4}-1)\text{(38)}$$

Dans ce cas le module est donc strictement supérieur à 1 et le schéma est instable.

On en déduit que la méthode d'Euler asymétrique est stable pour l'oscillateur harmonique si hω<2.