MENUCours D'introduction à l'analyse numérique

Il est fréquent de vouloir étudier un système dynamique conservatif sur une longue durée. Par exemple, la question de la stabilité du système solaire - dont il est nul besoin de rappeler l'importance - nécessite d'étudier sur des périodes séculaires un système quasi conservatif de 8 planètes et de leurs satellites. La solution analytique étant inaccessible, on met en place un traitement numérique du problème. Les algorithmes classiques d'ordre élevé (Runge-Kutta d'ordre 4) ont tendance à fournir une dynamique aux temps courts de bonne qualité mais ont le défaut de produire une dérive de l'énergie aux temps longs. Dans ce cas, on a recours à des méthodes numériques dites symplectiques, particulièrement adaptées aux systèmes conservatifs et supérieures aux méthodes classiques dans le sens où elles conduisent à une dérive de l'énergie faible aux temps longs.

Défaut des méthodes classiques

Une des raisons de l'inefficacité des méthodes classiques (Euler, Runge-Kutta) est leur caractère non réversible en temps : changer le pas temporel \(h\) en \(-h\) donne des équations différentes. Or, un système conservatif est réversible.

Voyons pourquoi la méthode d'Euler est intrinsèquement irréversible sur l'exemple d'une particule soumise à une force conservative :

La méthode d'Euler permet de ramener cette équation différentielle à un schéma itératif :

Changeons le sens d'écoulement du temps (\(h\to -h\) et \(n+1\to n-1\)). On obtient alors : \[\left\lbrace\begin{array}{ccc} \overrightarrow{v_n} & = & \overrightarrow{v}_{\!n-1} + \overrightarrow{a_n}\,h\\[3mm] \overrightarrow{r}_{n} & = & \overrightarrow{r}_{\!n-1} + \overrightarrow{v_n}\,h\\ \end{array}\right.\] ce qui ne correspond pas à ce que donne le schéma (1) à l'instant \(n\,h\). Le schéma d'Euler est donc irréversible alors que l'équation du mouvement ne l'est pas.

L'algorithme de Verlet à deux pas

Un des algorithmes symplectiques les plus simples à mettre en œuvre et largement utilisé en dynamique moléculaire est l'algorithme de Verlet. Il repose sur le développement de Taylor du vecteur position à l'ordre 3 aux instants \(t+h\) et \(t-h\) : \[\left\lbrace\begin{array}{ccc} \overrightarrow{r}(t+h) & = & \overrightarrow{r}(t)+h\overrightarrow{v}(t)+\dfrac{h^2}{2}\overrightarrow{a}(\overrightarrow{r}(t))+\dfrac{h^{3}}{3!}\dfrac{\text{d}^{3}\overrightarrow{r}}{\text{d}t^{3}}+\mathcal{O}(h^{4})\\[3mm] \overrightarrow{r}(t-h) & = & \overrightarrow{r}(t)-h\overrightarrow{v}(t)+\dfrac{h^2}{2}\overrightarrow{a}(\overrightarrow{r}(t))-\dfrac{h^{3}}{3!}\dfrac{\text{d}^{3}\overrightarrow{r}}{\text{d}t^{3}}+\mathcal{O}(h^{4}) \end{array}\right.\] En sommant ces deux équations, la vitesse disparaît et l'on obtient, aux erreurs d'ordre 4 près, le schéma itératif suivant : \[ \overrightarrow{r}_{\!n+1} = 2 \overrightarrow{r_{n}} -\overrightarrow{r}_{\!n-1}+ \overrightarrow{a_{n}}h^{2} \] Il est facile de constater que ce schéma est réversible en temps et donc particulièrement adapté à l'étude des systèmes conservatifs dont les forces ne dépendent que de la position. En revanche, cet algorithme présente deux défauts.

• Il s'agit d'un schéma à deux pas : l'itération ne peut démarrer que si l'on connaît \(\overrightarrow{r_{0}}\) et \(\overrightarrow{r}_{\!-1}\). Or, en général, les conditions initiales se résument par la donnée de \(\overrightarrow{r_{0}}\) et \(\overrightarrow{v}_{0}\). Il est d'usage alors d'initier l'itération à l'aide d'un développement de Taylor : \[\overrightarrow{r}_{\!\!-1}=\,\overrightarrow{r_0} - \overrightarrow{v_{0}}h+\frac{h^2}{2}\overrightarrow{a_0}\]
• Si l'on s'intéresse à l'évolution d'une grandeur faisant intervenir la vitesse (comme l'énergie cinétique), celle-ci se calcule à l'aide de la formule \[\overrightarrow{v_{n}}=\dfrac{\overrightarrow{r}_{\!n+1} - \overrightarrow{r}_{\!n-1}}{2h}\] On voit donc que le calcul de la vitesse consiste à soustraire deux nombres voisins ce qui produit des erreurs d'arrondi importants.

C'est pourquoi, on utilise en général une autre version algorithmique qui est mathématiquement équivalente à la version originale de Verlet mais qui présente l'intérêt d'être à un pas

L'algorithme de Verlet à un pas

L'algorithme de Verlet à un pas repose sur le schéma numérique suivant : \[\left\lbrace\begin{array}{ccc} \overrightarrow{r}_{\!n+1} & = & \overrightarrow{r_{n}} + h \overrightarrow{v_{n}}+\dfrac{h^2}{2}\overrightarrow{a_{n}}\\ \overrightarrow{v}_{\!n+1} & = & \overrightarrow{v_{n}} +\dfrac{h}{2}(\overrightarrow{a_{n}}+\overrightarrow{a}_{n+1})\\ \end{array}\right.\] On procédera donc ainsi :

Illustration avec le pendule simple

Considérons un pendule simple de masse \(m\), de longueur \(\ell\) dans un champ de pesanteur \(\overrightarrow{g}\). L'écart angulaire \(\theta\) par rapport à la position d'équilibre vérifie l'équation différentielle \[\ddot \theta=-\omega_{0}^{2}\sin \theta \qquad\text{avec}\qquad \omega_{0}^{2}=\sqrt{g/\ell}\] 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{ou}\qquad \left\lbrace\begin{array}{rcl} \dot\theta &=&\omega\\ \dot\omega &=&-\sin\theta\\ \end{array}\right.\]

Comparons les schémas numériques d'Euler et de Verlet (version à un pas) associés à cette équation différentielle. Posons \(\theta_{n}\) et \(\omega_{n}\) l'angle et la vitesse angulaire à l'instant \(t_{n}=nh\).

La méthode d'Euler donne \[\left\lbrace \begin{array}{c} \theta_{n+1}=\theta_{n}+\omega_{n}h\\ \omega_{n+1}=\omega_{n}-\sin\theta_{n}h\\ \end{array}\right.\] L'énergie mécanique du pendule simple est proportionnelle à \(E_{n}=\frac{1}{2}\omega_{n}^{2}-\cos{\theta_{n}}\). Le pendule simple étant un système conservatif, \(E_n\) doit rester constant. Or, on trouve \[E_{n+1} = E_{n}+h^{2}\left(\frac12\omega^{2}_{n}\cos\theta_{n}+\frac12\sin^{2}\theta_{n}\right)\] Le terme \(\left(1/2\omega^{2}_{n}\cos \theta_{n}+1/2\sin^{2} \theta_{n}\right)\) est en moyenne positif de sorte que l'on observe une dérive positive à long terme.

Quant à la méthode de Verlet à un pas, on a le schéma numérique suivant : \[\left\lbrace \begin{array}{c} \theta_{n+1}=\theta_{n}+\omega_{n}h-1/2\sin\theta_{n}h^{2}\\ \omega_{n+1}=\omega_{n}-1/2(\sin\theta_{n}+\sin\theta_{n+1})h \end{array}\right.\] À partir de cette rélation de récurrence on montre que \[E_{n+1}=E_{n}+\mathcal{O}(h^{3})\] En d'autres termes, à l'ordre 2, la méthode de Verlet respecte la conservation de l'énergie mécanique contrairement à la méthode d'Euler.

La simulation ci-dessous illustre bien le phénomène de dérive de la méthode d'Euler et la grande stabilité de celle de Verlet.