Intégrateur RK4 symplectique à pas fixe

Définition d'un intégrateur symplectique

On cherche à résoudre numériquement l'équation : $\dot{y} = f(y,t)$ . Un intégrateur numérique est dit symplectique s'il conserve la structure hamiltonienne des équations, c'est-à-dire les expressions des $\dot{q_i}$ et des $\dot{p_i}$ . Il conserve en particulier le Hamiltonien

(qui représente l'énergie du système) et donc l'intégrale de Jacobi. Un choix de conditions initiales du système place la particule étudiée sur une variété d'énergie précise dans l'espace des phases. L'intégrateur symplectique nous assure alors de rester sur celle-ci.

Définition du RK4 symplectique à pas fixe

Un tel intégrateur se construit en choisissant judicieusement les points où l'on évalue la fonction

entre deux itérations. Ces derniers sont au nombre de trois (différents du point de départ y(t)

), et sont notés $\xi_1, \xi_2$ et $\xi_3$ .
On passe alors du point y(t)

au point

via la relation :
$y(t+h) = y(t) + h \, b_1 f(t + \frac{b_1}{2} h, \xi_1) + h \, b_2 f(t + (b_1 + \frac{b_2}{2}) h, \xi_2) + h \, b_3 f(t + (b_1 + b_2 + \frac{b_3}{2}) h, \xi_3)$
où les $\xi_i$ sont donnés par les expressions suivantes :
$\xi_1 = y(t) + h \, \frac{b_1}{2} f(t + \frac{b_1}{2} h, \xi_1)$ ,
$\xi_2 = y(t) + h \, b_1 f(t + (b_1 + \frac{b_2}{2}) h, \xi_1) + h \, \frac{b_2}{2} f(t + (b_1 + \frac{b_2}{2}) h, \xi_2)$ ,
$\xi_3 = y(t) + h \, b_1 f(t + (b_1 + b_2 + \frac{b_3}{2}) h, \xi_1) + h \, b_2 f(t + (b_1 + b_2 + \frac{b_3}{2}) h, \xi_2) + h \, \frac{b_3}{2} f(t + (b_1 + b_2 + \frac{b_3}{2}) h, \xi_3)$
et où les b_i

valent (pour la méthode choisie; réf : "Numerical Hamiltonian Problems", J.M. Sanz-Serna et M.P. Calvo, p. 102) :
$b_1 = \frac{1}{3}(2 + 2^{1/3} + 2^{-1/3})$ ,
b_2 = 1 - 2 b_3

Les $\xi_i$ sont donnés implicitement; on parle alors de méthode implicite. Précisons qu'un intégrateur symplectique est forcément implicite (ibid.). La méthode est dite "DIRK" (Diagonally Implicit Runge-Kutta methods) : la relation définissant $\xi_1$ ne fait apparaître ni $\xi_2$ ni $\xi_3$ , et celle définissant $\xi_2$ ne fait pas intervenir $\xi_3$ . On peut donc d'abord déterminer $\xi_1$ , puis $\xi_2$ , et enfin $\xi_3$ .

Implémentation

Les coefficients $\xi_i$ ont été calculé par l'algorithme de Newton "classique". Pour le système étudié ici, la fonction

(l'équation du mouvement) est autonome, c'est-à-dire qu'elle ne dépend pas du temps. L'algorithme de Newton a donc d'abord été appliqué à la fonction
$\begin{array}{rcl} g_1 \, : \, \mathbb{R}^4 & \longrightarrow & \mathbb{R}^4 \\ x & \longmapsto & y(t) + h \, \frac{b_1}{2} f(x) - x \end{array}$
$\xi_1$ est alors la solution (sous réserve d'existence et d'unicité) de $g_1(\xi_1)=0$ .
On détermine ensuite $\xi_2$ comme zéro de la fonction
$\begin{array}{rcl} g_2 \, : \, \mathbb{R}^4 & \longrightarrow & \mathbb{R}^4 \\ x & \longmapsto & y(t) + h \, b_1 f(\xi_1) + h \, \frac{b_2}{2} f(x) - x \end{array}$
puis $\xi_3$ comme zéro de
$\begin{array}{rcl} g_3 \, : \, \mathbb{R}^4 & \longrightarrow & \mathbb{R}^4 \\ x & \longmapsto & y(t) + h \, b_1 f(\xi_1) + h \, b_2 f(\xi_2) + \frac{b_3}{2} f(x) - x \end{array}$
On peut alors calculer y(t+h)

. En réitérant ce processus, on peut alors déterminer la trajectoire complète de la particule.

Intégration numérique (2/3)