MENUCours D'introduction à l'analyse numérique

La méthode d'EulerMéthode inventée par le mathématicien Leonhard Euler en 1768. est une procédure numérique qui permet de résoudre de façon approximative des équations différentielles ordinaires du premier ordre avec condition initiale. Elle a le mérite d'être simple à comprendre et à programmer.

On cherche donc une solution approchée d'une équation ordinaire se mettant sous la forme \[\left\{\begin{array}{lccc} \dot{\mathbf{y}} & = & f(t,\mathbf{y}(t)), & 0\leq t \leq T \\ \mathbf{y}(0) & = & \mathbf{y}_{0} & \\ \end{array}\right.\] où \(\mathbf{y}(t)\) est un scalaire ou un vecteur et \(f(t,\mathbf{y}(t))\) une fonction suffisamment régulière pour que l'existence et l'unicité de la solution ne pose pas de problème.

Généralités

L'approche numérique

Les méthodes numériques employées pour résoudre les équations différentielles sont des méthodes approximatives basées sur la discrétisation de la variable \(t\) ainsi que sur l'utilisation de différences finies pour approcher les dérivées. Le problème se ramène alors à un calcul itératif, facile à automatiser à l'aide d'un programme informatique.

Pour effectuer ce calcul numérique, l'utilisateur doit disposer :

1. de la durée \(T\) de la simulation numérique ;
2. des conditions initiales et de la fonction \(f\) ;
3. du pas de discrétisation \(h\). L'intervalle \([0,T]\) est alors divisé en \(N\) subdivisions de même longueur \(h\). Finalement, la méthode numérique renvoie une liste \((y_{0},y_{1},,\ldots,y_{N})\) contenant les valeurs approchées de \(y(t_{n})\) pour les différents instants \(t_{n}=nh\).

On voit immédiatement que si l'on cherche le comportement de \(y(t)\) à long terme, le nombre de points à calculer peut devenir très important ce qui exige rapidité de calcul et de la mémoire. Insistons sur le fait que toutes les méthodes numériques ont leur point faible. Le bon choix est souvent un compromis entre simplicité de mise en œuvre, rapidité, stabilité et précision. Cette dernière est en effet limitée par les erreurs liées à l'algorithme (on parle d'erreurs de troncature) et celle liées à la machine (erreurs d'arrondi). En général, quand les unes diminuent les autres augmentent de telle sorte qu'il est impossible de minimiser, en même temps, et l'erreur d'arrondi et l'erreur de troncature : le meilleur choix est le fruit d'un compromis.

Travail préliminaire

Avant de résoudre numériquement une équation différentielle ordinaire, il est conseillé de procéder à un travail de simplification qui passe par la définition d'un nouveau système d'unités adapté au problème. En effet, pour minimiser les erreurs d'arrondi, il est bon de faire en sorte que les valeurs numériques utilisées soient de l'ordre de l'unité. Par exemple, pour traiter le problème de deux astres liés par la gravitation, il est préférable d'adopter le système d'unités astronomiques (le temps exprimé en années, les masses en masses solaires et les distances en rayon orbital terrestre) plutôt que le Système international (seconde, kilogramme et mètre). Or, changer d'unités, c'est finalement, définir de nouvelles grandeurs adimensionnées.

Prenons l'exemple de la chute libre d'une bille subissant une résistance proportionnelle à la vitesse régie par l'équation différentielle (relation fondamentale de la dynamique) \[m\frac{\mathrm{d} v}{\mathrm{d} t} = -\alpha v + mg \quad \text{avec}\quad v(0) = 0\] où \(m\) désigne la masse, \(\alpha\) le coefficient de frottement et \(g\) le champ de pesanteur. À partir de ces trois paramètres, on peut définir un temps caractéristique (temps de relaxation) \(\tau=m/\alpha\) et une vitesse caractéristique (vitesse limite) \(v_{\infty}=g\tau\). Ainsi, il est tentant d'adopter un nouveau système d'unités dans lequel :

1. le temps est mesuré en unité de \(\tau\), d'où la nouvelle variable temporelle \(t^{*}=t/\tau\) ;
2. la vitesse est mesurée en unité de \(v_{\infty}\) d'où la nouvelle variable \(v^{*}=v/v_{\infty}\) ;
3. la distance étant le produit d'une vitesse par un temps, est mesurée en unité de \(\ell\) avec \(\ell=v_{\infty}\tau\).

Avec ce nouveau choix d'unités et de variables adimensionnées, l'équation différentielle initiale prend une forme simple : \[\frac{\mathrm{d}v^{*}}{\mathrm{d}t^{*}} = 1- v^{*}\quad \text{avec}\quad v^{*}(0) = 0 \] Notez que cette équation différentielle est celle que nous aurions obtenue en faisant \(m=1\) et \(g=1\) et \(\alpha=1\). Finalement, on retiendra que lorsque l'on rend unitaire des paramètres physiques, cela équivaut à définir un nouveau jeu d'unités Il faut juste prendre garde à rendre unitaire des paramètres de dimensions indépendantes..

Enfin, ce qui va conditionner la valeur du pas \(h\) c'est l'ensemble des temps caractéristiques du phénomène que l'on modélise. Par exemple, dans le problème à \(N\) corps en orbite autour d'une étoile fixe, on peut définir \(N\) temps caractéristiques (les périodes orbitales par exemple). Si l'on veut décrire la dynamique du système de façon complète et précise il est impératif que, d'une part le pas d'intégration \(h\) soit plus petit que le plus petit des temps caractéristiques, d'autre part que la la durée \(T\) de la simulation soit supérieur au plus grand des temps caractéristiques. \[h\ll \tau_{\rm min} \quad\text{et}\quad T>\tau_{\rm max} \] Dans l'exemple de la chute libre, l'unique temps caractéristique vaut \(\tau=1\) dans le nouveau système d'unité. Ainsi, on prendra \(h\ll 1\) et \(T>1\). Par exemple, si l'on choisit \(h=0,01\) et \(T=10\), le nombre de points à calculer est \(N=10/0,01=1000\).

La méthode d'Euler

Algorithme d'Euler

La méthode présentée ici est dite méthode d'Euler explicite. Considérons l'équation différentielle ordinaire suivante \[ \left\{ \begin{array}{lclr} \dot{\mathbf{y}} & = & f(t,\mathbf{y}(t)), & 0\leq t \leq T \\ \mathbf{y}(0) & = & \mathbf{y}_{0} & \\ \end{array} \right. \] On peut intégrer cette équation comme suit : \[ \mathbf{y}(t_{n+1})=\mathbf{y}(t_{n})+\int_{t_{n}}^{t_{n+1}}f(t,\mathbf{y}(t))\;\mathrm{d}t \] La méthode d'Euler consiste à approcher l'intégrale par la méthode des rectangles à gauche : \[\int_{t_{n}}^{t_{n+1}}f(t,\mathbf{y})\;\mathrm{d}t \simeq h\times f(t_{n},\mathbf{y}(t_{n}))\] D'où le schéma itératif suivant

Exemple d'une équa-diff d'ordre un

Appliquons cette méthode au problème de la chute libre décrite au §1.2. L'équation différentielle s'écrit \[\frac{\mathrm{d} v^{*}}{\mathrm{d} t^{*}} = 1- v^{*}\quad \text{avec}\quad v^{*}(0) = 0\] et le schéma itératif (1) donne : \[\left\{\begin{array}{lcl} v^{*}_{n+1} & = & v^{*}_{n} + h(1-v^{*}_{n}) \\[1mm] v^{*}_{0} & = & 0 \\ \end{array}\right.\]

On peut voir ici que la méthode d'Euler est une approximation au premier ordre en \(h\). En effet, calculons le premier terme \(v^{*}_{1}\) et comparons le au terme exacte \(v^{*}(h)\), sachant que la solution exacte est \(v^{*}(t^{*})=1-e^{-t^{*}}\). On a \[v^{*}_{1} = v^{*}_{0} + h(1-v^{*}_{0})= h \] à comparer avec \[v^{*}(h)=1-e^{-h}\] Il apparaît immédiatement que \(v^{*}_{1}\) et \(v(h)\) concordent si l'on remplace \(e^{-h}\) par son développement limité au premier ordre, à savoir \(e^{-h}=1-h+{\cal}O (h^{2})\).

La figure 1 montre la solution numérique pour différents pas ainsi que la solution exacte. On constate que plus le pas \(h\) est petit meilleure est l'adéquation avec la solution exacte. On constate qu'avec un pas \(h=0,2\), l'erreur est visible à l'œil nu (de l'ordre de quelques %). Ça n'est que lorsque le pas est très petit devant un, que l'erreur est insignifiante.

Choisissons maintenant \(h=2\), autrement dit un pas de discrétisation plus grand que le temps caractéristique \(\tau\) (qui vaut 1 puisque \(\tau\) est l'unité de temps). On voit alors apparaître un phénomène d'instabilité numérique (Fig. 2). C'est l'un des inconvénients de la méthode d'Euler : la stabilité n'est pas toujours garantie. C'est pourquoi il est recommandé de tester différents pas pour vérifier que le résultat est robuste. Dans ce cas, cela signifie généralement que l'on est loin des conditions d'instabilité.

Exemple d'une équa-diff d'ordre deux

La plupart du temps, l'évolution d'un système physique obéit à une équation différentielle qui n'a pas forcément le bon goût d'être une équation différentielle scalaire du premier ordre. On rencontre plus souvent des systèmes d'équations différentielles ou une équation différentielle scalaire d'ordre supérieure à un. En mécanique, il arrive couramment que la dynamique d'un système soit décrite par une équation différentielle d'ordre deux.

Par exemple, considérons le mouvement d'un pendule simple de masse \(m\) et de longueur \(\ell\) qu'on lâche dans le champ de pesanteur \(g\) après l'avoir écarté de l'angle \(\theta_0\). L'équation du mouvement s'écrit \[\ddot \theta=-{\omega_{0}}^{2}\sin \theta \qquad\text{avec}\qquad \left\{\begin{array}{lcl} {\omega_{0}}^{2}&=&\sqrt{g/\ell}\\ \theta(0) &=&\theta_0\\ \dot \theta(0) &=&0\\ \end{array}\right. \] Tout d'abord, changeons l'unité de temps de façon à simplifier l'équation différentielle. Posons \(\tau=\sqrt{\ell/g}\) notre nouvelle unité de temps ce qui revient à redéfinir le temps par \(t_{*}=t/\tau\). Il est alors facile de montrer qu'après cette redéfinition du temps, l'équation différentielle se simplifie \[\ddot\theta=-\sin\theta \qquad\text{avec}\qquad \left\{\begin{array}{rcl} \theta(0) &=&\theta_0\\ \dot \theta(0) &=&0\\ \end{array}\right. \]

Cette équation peut se mettre sous la forme d'une équation différentielle d'ordre un, à condition de l'écrire sous forme vectorielle. En effet, définissons le vecteur de dimension 2 dont la première composante est l'angle \(\theta\) et la seconde la vitesse angulaire \(\omega=\dot \theta\) : \[ \mathbf{y}= \left(\begin{array}{c}\theta\\\omega\end{array}\right) \quad\text{avec}\quad \mathbf{y}(0)= \left(\begin{array}{c}\theta_0\\ 0\end{array}\right) \] On peut alors écrire l'équation du pendule simple sous la forme \[ \frac{\mathrm{d}}{\mathrm{d}t} \left(\begin{array}{c}\theta\\ \omega\end{array}\right) =\left(\begin{array}{c}\omega \\-\sin \theta\end{array}\right) =f(t,\mathbf{y}(t)) \] Posons \(\theta_{n}\) et \(\omega_{n}\) l'angle et la vitesse angulaire à l'instant \(t_{n}=nh\). La méthode d'Euler donne

On peut constater, sur le résultat de la simulation, un autre défaut de la méthode d'Euler : bien que le système étudié soit conservatif, le schéma d'Euler ne respecte pas la conservation de l'énergie. En effet, en vertu des lois de la mécanique, l'amplitude des oscillations doit rester constante, contrairement à ce que l'on observe. Plus précisément, l'énergie mécanique du pendule simple à l'instant \(t_n\) est proportionnelle à \(E_{n}=\frac{1}{2}\omega_{n}^{2}-\cos{\theta_{n}}\). Or, à partir de la récurrence (2), on obtient \[E_{n+1} = E_{n}+h^{2}\left(\frac12\omega^{2}_{n}\cos\theta_{n}+\frac12\sin^{2}\theta_{n}\right)\] Par conséquent, puisque le terme \(\left(1/2\omega^{2}_{n}\cos \theta_{n}+1/2\sin^{2} \theta_{n}\right)\) est positif en moyenne, il en découle une dérive positive de l'énergie quel que soit le pas. Autrement dit, le comportement à long terme est mal décrit. La simulation suivante illustre cette dérive.

En pratique, lorsque l'on met en place une résolution numérique d'un système conservatif, on fait plutôt appel à l'algorithme de Verlet[1] particulièrement adapté dans le sens où elle conduit à une dérive de l'énergie faible aux temps longs.

Conclusion

Bien que la méthode d’Euler présente l'inconvénient de propager facilement les erreurs, et même de les amplifier, elle n'est pas dénuée d'intérêt : simple à programmer, elle produit rapidement une solution approximative qui, si le pas est bien choisi et la durée d'observation raisonnable, permet d'avoir une première approche du phénomène étudié.