F

MENUCours Outils et Méthodes pour la Physique

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.

Algorithmes de Verlet

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 :

\begin{equation} \left\{\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.\label{eq:C5eqdumouvement} \end{equation}

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

\begin{equation} \left\{\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. \label{eq:EulerIrreversible} \end{equation}

Changeons le sens d'écoulement du temps ($h\to -h$ et $n+1\to n-1$). On obtient alors : \[\left\{\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 \eqref{eq:EulerIrreversible} à 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\{\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 : \begin{equation} \overrightarrow{r}_{\!n+1} = 2 \overrightarrow{r_{n}} -\overrightarrow{r}_{\!n-1}+ \overrightarrow{a_{n}}h^{2} \label{eq:C4verlet} \end{equation} 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 :

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\{\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

  1. Initialisation du pas $h$, de la durée $T$.
  2. Initialisation des conditions initiales : $t=0$, $\overrightarrow{r}=\overrightarrow{r_{0}}$ et $\overrightarrow{v}=\overrightarrow{v_{0}}$.
  3. Définition de la fonction $\overrightarrow{a}(\overrightarrow{r})$.
  4. Tant que $t\leq T$ faire :
    1. Calcul de $\overrightarrow{a_0}=\overrightarrow{a}(\overrightarrow{r})$.
    2. Nouvelle position : $\overrightarrow{r}=\overrightarrow{r}+h\overrightarrow{v}+\frac{h^2}{2}\overrightarrow{a_0}$.
    3. Calcul de $\overrightarrow{a_1}=\overrightarrow{a}(\overrightarrow{r})$.
    4. Nouvelle vitesse : $\overrightarrow{v}=\overrightarrow{v}+\frac{h}{2}(\overrightarrow{a_0}+\overrightarrow{a_1})$
    5. $t=t+h$.
    6. Enregistrement des données

Exemples

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\{\begin{array}{rcl} \dot\theta &=&\omega\\ \dot\omega &=&-\sin\theta\\ \end{array}\right.\]

Remarque

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\{ \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\{ \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.

Simulation

Built with Processing

Evolution d'un pendule simple lâché sans vitesse avec un angle de 45°. La simulation dure 1 minute.

Your browser does not support the canvas element.

Le problème à deux corps

La méthode de Runge-Kutta d'ordre 4 si utilisée présente également ce défaut de ne pas conserver l'énergie mécanique [1]. On peut le voir dans le problème à deux corps suivant.

Supposons deux astres de masse $m_{1}=1~m_{\star}$ et $m_{2}=0.1~m_{\star}$ en interaction gravitationnelle. Initialement séparés de $2,5\;\mathrm{ua}$, ils sont lancés avec des vitesses de sens opposé : $v_{01}=-0,3\;\mathrm{u.a/an}$ et $v_{02}=3\;\mathrm{u.a/an}$. Notez que la quantité de mouvement du système étant nulle, le centre de gravité reste fixe.

Sur le plan théorique, on montre que les deux astres doivent décrire des ellipses de foyer le centre de gravité du système double. On peut calculer que l'astre le plus léger doit décrire une ellipse de grand-axe $a=3,64\;\mathrm{ua}$ et de période $T=2,34\;\mathrm{ans}$. L'équation à laquelle obéissent les vecteurs positions des deux astres s'écrit dans le système d'unités astronomiques : \[\dfrac{\mathrm{d}^{2}\overrightarrow{r_{i}}}{\mathrm{d}t^{2}}=-4\pi^{2}\frac{m_{j}(\overrightarrow{r_{i}}- \overrightarrow{r_{j}})}{r_{ij}^{3}}\] La figure ci-dessous montre les trajectoires que l'on obtient à l'aide des méthodes de Runge Kutta et de Verlet. On note que la méthode de Runge Kutta produit une dérive énergétique qui se traduit par une trajectoire en spirale : aux temps longs, l'astre de masse $m_{2}$ va finir par être éjecté du système ce qui contredit complètement le résultat exacte qui prévoit une trajectoire périodique stable.

trajectoire obtenue par la méthode de Verlet trajectoire obtenue par la méthode de Runge-Kutta d'ordre 4
Simulations numériques de la trajectoires de deux astres liés par gravitation avec un pas $h=0,01\;\mathrm{an}$. Les points correspondent à une photographie toutes les 0,08 ans.

Vous aimez ?

Pour en savoir plus...

  1. Roussel J. Méthodes de Runge-Kutta. FEMTO, la physique enseignée http://femto-physique.fr/omp/runge_kutta.php