Intégrateur symplectique

Ces intégrateurs ont pour but de résoudre les équations d'Hamilton, qui s'écrivent :

${\displaystyle {\dot {q}}_{i}={\frac {\partial {\mathcal {H}}}{\partial p_{i}}}\quad {\mbox{et}}\quad {\dot {p}}_{i}=-{\frac {\partial {\mathcal {H}}}{\partial q_{i}}}.}$

où ${\displaystyle {\boldsymbol {q}}=(q_{1},\dots ,q_{N})\in \mathbb {R} ^{N}}$ sont les coordonnées généralisées, ${\displaystyle {\boldsymbol {p}}=(p_{1},\dots ,p_{N})\in \mathbb {R} ^{N}}$ les moments généralisés conjugués de ces coordonnées et ${\displaystyle {\mathcal {H}}}$ est l'Hamiltonien du système. Le point au-dessus de ces variables dénote la dérivée temporelle de celles-ci. Toutes les variables ${\displaystyle ({\boldsymbol {p}},{\boldsymbol {q}})}$ respectant ces équations sont dites canoniques.

Par construction, l'évolution temporelle de ces variables conserve la 2-forme symplectique.

${\displaystyle \omega =\sum _{i=1}^{N}dp_{i}\wedge dq_{i}=\sum _{i=1}^{N}dp'_{i}\wedge dq'_{i}}$

où ${\displaystyle ({\boldsymbol {q}}',{\boldsymbol {p}}')}$ sont les variables à un instant ultérieur et l'opérateur ${\displaystyle \wedge }$ est un produit extérieur. Cette transformation est donc un symplectomorphisme. Un schéma d'intégration est symplectique lorsqu'il préserve cette 2-forme à chaque pas de temps.

D'après le Théorème de Liouville, le volume dans l'espace des phases le long du flot de ${\displaystyle {\mathcal {H}}}$ est conservé. Un intégrateur symplectique va donc lui aussi conserver ce volume, c'est une propriété très intéressante lorsqu'on veut faire des intégrations sur des temps longs par rapport à l'échelle de temps du système étudiée. Cela induit qu'un nuage de conditions initiales ne peut ni se contracter ni se dilater^[2], contrairement aux schémas classiques qui peuvent introduire des attracteurs inexistants dans le système hamiltonien. C'est particulièrement important en mécanique céleste pour maintenir la stabilité des éléments orbitaux sur des temps de plusieurs millions d'années.

La majorité des intégrateurs usuels, comme la méthode d'Euler ou les méthodes de Runge-Kutta, ne sont pas des intégrateurs symplectiques.

Démonstration pour la méthode d'Euler

On cherche à montrer que la méthode d'Euler n'est pas symplectique à partir du cas du pendule simple. L'hamiltonien du système se note :

${\displaystyle {\cal {H={\frac {p^{2}}{2}}+\cos q}}}$

Ainsi, d'après les équations canoniques d'Hamliton on obtient le système suivant.

${\displaystyle \left\{{\begin{array}{rcl}{\dot {p}}&=&\sin q\\{\dot {q}}&=&p\end{array}}\right.}$

En notant ${\displaystyle (q',p')}$ les variables du mouvements au temps ${\displaystyle t+\tau }$. La méthode d'Euler s'écrit :

${\displaystyle \left\{{\begin{array}{rcl}p'&=&p+\tau \sin q\\q'&=&q+\tau p\end{array}}\right.}$

On montre donc rapidement.

${\displaystyle dp'\wedge dq'=d(p+\tau \sin q)\wedge d(q+\tau p)=(1-\tau ^{2}\cos q)dp\wedge dq\neq dp\wedge dq}$

La 2-forme n'est pas conservée par la méthode d'Euler, ce schéma d'intégration n'est donc pas symplectique.

Souvent, H = A + ε B, où A est intégrable et B est une perturbation intégrable souvent aussi, et ε un réel très petit. Appelons L l'opérateur de Liouville de A et M l'opérateur de Liouville de B :

Alors le problème est de calculer exp{L + ε M} t qui hélas est différent de exp{L t} .(exp Mt)^ε.

Le cas classique en mécanique céleste est la perturbation de Saturne par Jupiter. Mais on peut aussi bien tester la méthode sur une particule dans un puits de potentiel.

Évidemment, il y a deux possibilités : la formule de Trotter ou la formule de Campbell-Hausdorff.

Ou bien des formes raffinées de combinaison des deux adéquates.

L'idée forte est la suivante : t est petit ; faire une théorie au n-ième ordre, conduit à une erreur O(tⁿ ε) ou plus exactement O(tⁿ ε +t² ε²) : on a intérêt à pousser la méthode jusqu'à l'ordre n tel que :

t^(n-2) = ε.

Dans certains cas, dits symétriques, on peut l'améliorer en t^(n-4) = ε.

On peut commencer par tester la méthode sur l'oscillateur harmonique, qui est le test usuel.

On continuera avec le pendule simple.

Cette fois, en coordonnées réduites, A = p²/2 et B = 1 - cos q . Le système est intégrable exactement, via les fonctions de Jacobi, mais nous préférons prendre A = p²/2 + ε q²/2 et B = 1 - cos q - q²/2

L'opérateur A est donc celui d'un oscillateur harmonique, et B joue le rôle de perturbation, si les oscillations sont pendulaires. Dans le cas de tournoiement, le problème de la "séparatrice" ne peut se régler simplement, car l'intégrateur symplectique ne conserve pas rigoureusement l'énergie. Il vaut mieux alors se tourner vers la solution du pendule simple discret.