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 : \[ \left\lbrace\begin{array}{ccl} \dfrac{\mathrm{d}\overrightarrow{v}}{\mathrm{d} t} &= & \dfrac{\overrightarrow{f}(\overrightarrow{r})}{m}=\overrightarrow{a}(\overrightarrow{r}) \\[2mm] \dfrac{\mathrm{d}\overrightarrow{r}}{\mathrm{d} t} & = & \overrightarrow{v} \\ \end{array} \right. \] La méthode d'Euler permet de ramener cette équation différentielle à un schéma itératif : \begin{equation} \left\lbrace\begin{array}{ccc} \overrightarrow{v}_{\!n+1} & = & \overrightarrow{v}_{\!n} + \overrightarrow{a}_{\!n}\,h\\[2mm] \overrightarrow{r}_{\!n+1} & = & \overrightarrow{r}_{\!n} + \overrightarrow{v}_{\!n}h\\ \end{array} \right. \end{equation} 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} - h\,\overrightarrow{v_{0}} +\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 :
Algorithme de Verlet à un pas
- Initialisation du pas \(h\), de la durée \(T\).
- Initialisation des conditions initiales : \(t=0\), \(\overrightarrow{r}=\overrightarrow{r_{0}}\) et \(\overrightarrow{v}=\overrightarrow{v_{0}}\).
- Définition de la fonction \(\overrightarrow{a}(\overrightarrow{r})\).
- Tant que \(t\leq T\) faire :
- Calcul de \(\overrightarrow{a_0}=\overrightarrow{a}(\overrightarrow{r})\).
- Nouvelle position : \(\overrightarrow{r} = \overrightarrow{r} + h\overrightarrow{v}+ \frac{h^2}{2}\overrightarrow{a_0}\).
- Calcul de \(\overrightarrow{a_1} = \overrightarrow{a}(\overrightarrow{r})\).
- Nouvelle vitesse : \(\overrightarrow{v} = \overrightarrow{v} + \frac{h}{2}(\overrightarrow{a_0} + \overrightarrow{a_1})\)
- \(t=t+h\).
- Enregistrement des données
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. \]
Notez que dans ce nouveau système d'unités, tout se passe comme si \(g=1\) et \(\ell=1\).
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.
Pour en savoir plus...
- Méthodes de Runge-Kutta[en ligne], 2015. Disponible sur femto-physique.fr