Sommaire
1. Enseignement Méthodologique: méthodes numériques
1. Problème à trois corps restreint 2007
1. Présentation du problème
1. Introduction
2. Obtention d'une formulation hamiltonienne adimensionnée et autonome
3. Formulation hamiltonienne du problème
4. Intégrale de Jacobi
2. Courbes de vitesse nulle
1. Caracterisation et détermination des courbes de vitesse nulle
2. Détermination numérique des courbes de vitesse nulles
3. Intégration numérique
1. Résumé
2. Intégrateur RK4 symplectique à pas fixe
3. RK4 classique à pas variable ou RK4 symplectique à pas fixe ?
4. Sections de Poincaré
1. Sections de Poincaré
5. Détermination du zéro d'une fonction
1. Méthodes de la dichotomie et de Newton-Raphson
2. Méthode de Newton multidimensionnelle
chapitre 1

Enseignement Méthodologique: méthodes numériques

chapitre 1, souschapitre 1

Problème à trois corps restreint 2007

chapitre 1, souschapitre 1, section 1

Présentation du problème

Introduction

Problème physique

On s'intéresse au problème à trois corps restreint, plus précisément au mouvement d'une particule test de masse m dans le champ de gravité de deux masses m_1 et m_2, m_2 < m_1, en rotation circulaire uniforme autour de leur centre de gravité commun. On suppose de plus que la particule test ne perturbe pas les deux masses, et qu'elle se meut dans le même plan qu'elles.

Modélisation numérique

La modélisation numérique a été faite en langage JAVA dans l'environnement de développement Simulab. Une applet a été réalisée et est disponible, ainsi que son mode d'emploi, sur ce site .

Points abordés

Dans ce chapitre, nous allons donner les parties de cours correspondant à la position et l'étude du problème, mais aussi à certains calculs effectués en arrière-plan lors de l'affichage de l'applet, tels que la recherche des points de Lagrange, les courbes de vitesses nulles, l'intégration numérique, et la détermination de zéro d'une fonction.



chapitre 1, souschapitre 1, section 1, page 1

Présentation du problème   (1/3)

Obtention d'une formulation hamiltonienne adimensionnée et autonome

Nous résumons ici les étapes nécessaires à l'obtention de la formulation hamiltonienne que nous utilisons par la suite. Le détail des calculs est donné ici. On écrit tout d'abord le Lagrangien associé au mouvement de la particule dans le champ gravitationnel créé par les deux masses m_1 et m_2. Afin de simplifier les expressions obtenues, on procède à un "dédimensionnement" des équations ; les constantes du problème sont alors incluses dans les variables et celles-ci deviennent dimensionnées. Bien que la formulation lagrangienne permette de décrire le système, la formulation hamiltonienne, plus adaptée aux problèmes de mécanique céleste, va être privilégiée ici. On dérive ainsi du Lagrangien le Hamiltonien associé au problème. Cependant, celui-ci dépend du temps. On le rend autonome (sans dépendance explicite du temps) via une transformation canonique (changement de variables qui préserve la structure hamiltonienne). chapitre 1, souschapitre 1, section 1, page 2

Présentation du problème   (2/3)

Formulation hamiltonienne du problème

Le Hamiltonien H finalement obtenu s'écrit :
H = \frac{1}{2} (P_1^2 + P_2^2) + P_1 Q_2 - P_2 Q_1 - \left(\frac{1- \mu}{R_1} + \frac{\mu}{R_2}\right)
avec comme équations associées :
\left\{\begin{array}{ll}\dot{Q_1} = \frac{\partial H}{\partial P_1} & \dot{P_1} = - \frac{\partial H}{\partial Q_1} \\\dot{Q_2} = \frac{\partial H}{\partial P_2} & \dot{P_2} = - \frac{\partial H}{\partial Q_2} \end{array}\right.
Les variables Q_1 et Q_2 se rapportent à la position de la particule, P_1 et P_2 sont leurs variables dites "conjuguées" et sont liées à la vitesse de la particule. Le Hamiltonien H représente en fait l'énergie du système. Théoriquement, il doit donc être conservé au cours du mouvement. On dit alors qu'il constitue une intégrale première du mouvement. Explicitement, les équations du mouvement sont donc :
\left\{\begin{array}{l}\dot{Q_1} = P_1 + Q_2 \\\dot{Q_2} = P_2 - Q_1 \\ \dot{P_1} = P_2 +(1-\mu) \cdot \frac{\partial}{\partial Q_1}(1/R_1) + \mu \frac{\partial}{\partial Q_2} (1/R_2) \\ \dot{P_2} = -P_1 +(1-\mu) \cdot \frac{\partial}{\partial Q_1}(1/R_1) + \mu \frac{\partial}{\partial Q_2} (1/R_2)\end{array}\right.

\left\{\begin{array}{l}R_1 = \sqrt{ (Q_1 + \mu )^2 +Q_2^2 } \\ R_2 = \sqrt{ (Q_1 + \mu - 1)^2 +Q_2^2 }\end{array}\right. chapitre 1, souschapitre 1, section 1, page 3

Présentation du problème   (3/3)

Intégrale de Jacobi

L'intégrale de Jacobi J est définie par J = -2 H, où H est le Hamiltonien du système. C'est donc une intégrale première du mouvement. Cela se traduit dans l'espace des phases (Q_1, Q_2, P_1, P_2) par le fait que le système doive rester sur la variété (de dimension 3) J = cste au cours du temps. Cette contrainte sur l'évolution est particulièrement importante lorsque l'on s'intéresse à l'évolution du système sur le long terme, par exemple lorsque l'on étudie des sections de Poincaré. L'intégrateur Runge-Kutta d'ordre 4 à pas variable utilisé ici ne conserve pas cette quantité. Sur l'applet, on observe effectivement une dérive constante de cette valeur. L'autre intégrateur disponible est de type Runge-Kutta d'ordre 4, à pas fixe, symplectique, qui lui conserve l'intégrale de Jacobi. On observe effectivement une plus grande stabilité de cette dernière sur les simulations de l'applet. Si l'on cherche à tracer des sections de Poincaré du système, il faudra donc utiliser l'intégrateur symplectique. chapitre 1, souschapitre 1, section 2, page 1

Courbes de vitesse nulle   (1/2)

Caracterisation et détermination des courbes de vitesse nulle

L'expression de l'Hamiltonien peut se réécrire sous une forme différente, en réintroduisant les vitesses relatives dans le référentiel en rotation accent(Q_1;.) et accent(Q_2;.), on obtient alors :
H=(1/2)*( (P_1+Q_2)^2 + (P_2-Q_1)^2 ) - (1/2)*(Q_1^2+Q_2^2)-((1-mu)/R_1 + mu/R_2 )
Puisque accent(Q_1;.) = P_1+Q_2et accent(Q_2;.)=P_2-Q_1, et que la courbe de vitesse nulle s'écrit
V^2=accent(Q_1;.)^2+accent(Q_2;.)^2,
nous obtenons :
V^2=accent(Q_1;.)^2+accent(Q_2;.)^2=(Q_1^2+Q_2^2)+2*((1-mu)/R_1 + mu/R_2) +2*H
Or, nous avons l'intégrale de Jacobi (lien) qui vaut J=-2H. L'équation de la courbe de vitesse nulle est donc la suivante :
V^2=((Q_1)^2+(Q_2)^2)+2*( (1-mu)/R_1 + mu/R_2) - J
Cette expression devant être toujours positive, on voit que pour une valeur fixée de J, le plan (Q1,Q2) sera partagé en plusieurs régions:

Bien entendu, ces courbes ne sont visible que dans le repère tournant.

Pour une même valeur du paramètre de masse (ici mu=0.1), il y a cinq types de comportements possibles :

chapitre 1, souschapitre 1, section 2, page 2

Courbes de vitesse nulle   (2/2)

Détermination numérique des courbes de vitesse nulles

Introduction

Les courbes de vitesse nulle dans le cadre de notre étude du problème à trois corps restreint sont des courbes fermées du plan ((Q_1*,Q_2)). Leur comportement topologique diffère suivant la valeur de J (voir chapitre précédent). Une fois ce comportement identifié, leur détermination numérique va donc se faire en trois étapes. On détermine tout d'abord un point de chacune des courbes que l'on considère, puis on suit le point trouvé le long de la courbe par continuation tout en définissant un paramètre de contrôle permettant l'arrêt de cette continuation (les courbes étant fermées).

Evaluation du comportement topologique des courbes de vitesse nulle suivant la valeur de J

Comme il est décrit au chapitre précèdent, il y a 5 comportement différents pour les courbes de vitesse nulle suivant la valeur de l'intégrale de Jacobi. Il convient donc de déterminer ces différents cas pour une détermination numérique appropriée des courbes.

Pour des conditions initiales données, et donc une intégrale de Jacobi J donnée, on rappelle que la courbe de vitesse nulle est définie par:
f(Q_1*,Q_2)=Q_1^2+Q_2^2+2*( (1-mu)/R_1 + mu/R_2) - J
avec système(R_1=((Q_1+mu)^2+(Q_2^2))^(1/2);R_1=((Q_1+mu-1)^2+(Q_2^2))^(1/2))
Les points de Lagrange L_4 et L_5 sont les premiers points où les courbes de vitesse nulle apparaissent. En se plaçant en L_4, on peut donc déterminer le seuil d'apparition des courbes. L_4 étant de coordonnées Q_1=(1/2-mu) et Q_2=sqrt(3)/2, on a R_1 et R_2 égaux à 1 et :
f(L_4)=3-mu+mu^2-J
Les courbes de vitesse nulle existent donc pour:
J>3-mu+mu^2
Les changement de topologies se font lorsque les différentes branches se touchent aux points de Lagrange. On peut donc déduire le nombre de branches à calculer en évaluant f en L_1, L_2 et L_3. On utilise bien sûr les valeurs numériques de Q_1 calculées lors de la détermination des points fixes, Q_2 étant nul pour ces points. Ainsi :

Recherche d'un point de départ pour la continuation d'une courbe

Afin de commencer la continuation de la courbe de vitesse nulle considérée, il nous faut d'abord trouver un point de cette courbe.
Pour cela on utilise les deux méthodes de détermination de zéro : dichotomie et Newton-Raphson dans cet ordre pour les raisons exposées dans les chapitres leur correspondant. En effet en prenant un point à l'intérieur de la courbe (où f<0) et un autre à l'extérieur (où f>0), il existe un point dans cet intervalle tel que f=0 et qui appartient donc à la courbe de vitesse nulle.
Pour la méthode par dichotomie on choisit les intervalles de départ selon les cas topologiques comme suit:
On démarre ensuite la méthode de Newton-Raphson depuis l'une des bornes de l'intervalle trouvé par dichotomie.

Continuation des courbes

Maintenant que l'on a trouvé un point de la courbe, on peut suivre le point trouvé le long de la courbe par "continuation". Pour cela écrivons un élément d'abscisse curviligne ds:
(ds)^2=(dQ_1)^2+(dQ_2)^2
or dQ_1 et dQ_2 sont reliés via la différentielle de f :
(partf/partQ_1)*dQ1+(partf/partQ_2)*dQ2=0
On déduit de ces deux équations le système suivant:
système(dQ_1/ds=-(partf/partQ_2)*(1/(sqrt(((partf/partQ_2))^2+((partf/partQ_1))^2)));dQ_2/ds=(partf/partQ_1)*(1/sqrt((((partf/partQ_2))^2+((partf/partQ_1))^2))))
On a ainsi un système de deux équations différentielles à deux inconnues donnant l'évolution de Q_1 et Q_2 en fonction de s pouvant être résolu avec un intégrateur numérique(lien). On a choisi dans le calcul pour l'applet un Runge-Kutta d'ordre 4 à pas variable.

Remarque :


Arrêt de la continuation

Le problème auquel on est confronté lors de la continuation des courbes de vitesse nulle est lié au fait que ces courbes sont fermé et que l'intégrateur calculant les points de la courbe ne peut testé intrinsèquement si il est revenu à son point de départ ou non. Il nous faut donc trouver un moyen de tester ce paramètre afin d'éviter au programme de calculer pour rien, et/ou de passer à une étape suivante.

On effectue donc un changement d'origine du repère ((Q_1*,Q_2))en le centrant sur un point à l'intérieur de la courbe, choisi astuscieusement en fonction des cas présentés plus haut. Pour différencier les axes de ce nouveau repère, on appellera x l'axe correspondant à Q_1 et y celui de Q_2. Ainsi, lorsque l'ordonnée est négative (y<0) et l'abscisse change de signe (i.e. x<0 devient x>0), un angle de 2*(pi) par rapport au premier point de la courbe est effectué. Autrement dit, la courbe se referme, on obtient une condition permettant de stopper la continuation.

chapitre 1, souschapitre 1, section 3, page 1

Intégration numérique   (1/3)

Résumé

Un cours de présentation de la méthode de Runge-Kutta (RK) "classique" est exposée ici Nous allons maintenant présenter la méthode de Runge-Kutta d'ordre 4 dite symplectique à pas fixe (RK4 symp). chapitre 1, souschapitre 1, section 3, page 2

Intégration numérique   (2/3)

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 H (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 f 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 y(t+h) 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,
b_3=b_1
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 f (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_1est 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.

chapitre 1, souschapitre 1, section 3, page 3

Intégration numérique   (3/3)

RK4 classique à pas variable ou RK4 symplectique à pas fixe ?

Les deux méthodes d'intégration ont chacune leurs avantages et inconvénients.
L'intégrateur RK4 à pas variable est plutôt utilisé lorsque l'on veut une grande précision sur des temps courts d'intégration. La dérive de l'intégrale de Jacobi est alors faible.
L'intégrateur RK4 symplectique est par contre à préférer lorsque l'on souhaite une précision constante sur le long terme (grand nombre d'itérations), et la conservation, en moyenne, de l'intégrale de Jacobi (application : sections de Poincaré). L'inconvénient d'un pas fixe se manifeste près des deux masses. Dans ces zones, si le pas d'intégration est trop grand, la méthode divergera. Il faut donc fixer à l'avance un pas d'intégration adapté, ce qui n'est pas toujours évident. De plus, un pas très petit consommera des ressources informatiques inutiles dans les zones où la trajectoire est peu courbée. chapitre 1, souschapitre 1, section 4, page 1

Sections de Poincaré   (1/1)

Sections de Poincaré

La section de Poincaré est souvent utilisé pour caractériser le chaos, et étudier la dynamique d'un système. Pour définir cette section de Poincaré dans le cadre du problème à trois corps restreint, rappelons nous qu'on s'intéresse à la particule test, dans le repère tournant. Il est alors intéressant de reporter sur un diagramme les valeurs (q1,p1) de cette particule à chaque fois qu'elle coupe l'axe q2=0. On dit alors que la section de Poincaré est à q2=0.

Plus généralement, une section de Poincaré se dessine dans un plan ayant pour abscisse la coordonnée généralisé q_{i} et pour ordonnée la coordonnée généralisé p_{i} (ou, quand le problème permet une formulation hamiltonienne en coordonnées action-angle : (angle-action)), p_{j} (j\nei) étant calculés pour que J=cste (on parle alors de variété d'énergie sur laquelle on se place, ie. une section de Poincaré est un ensemble de trajectoires à énergie constante), et q_{j} (j\nei) étant libre. La section est alors généralement une surface de section (ie. équation de plusieurs variables).
Dans notre cas, il est intéressant de tracer une section de Poincaré à q2=0.

Par ailleurs, on définit une trajectoire de Poincaré, les points de ce diagramme issus de même conditions initiales q10, q20, p10, p20. On fixe généralement q20 une fois pour toute, et on calcule p20 en fonction de q10, p10 et J (qui doit rester constant sur toutes les trajectoires). Plus le nombre de points d'une trajectoire de Poincaré est élevé, plus la section de Poincaré est claire à l'interprétation.

On obtient alors des courbes fermées (si le nombre de trajectoires de Poincaré est suffisant) correspondant à l'ensemble des valeurs des couples (q1,p1) prises par la particule test quand elle suit une orbite non-chaotique, et des zones chaotiques caractérisées par une densité de points non ordonnés élevée, et prises entre deux courbes fermées.
Les zones chaotiques entourent généralement un ensemble de courbes fermées, appelé île. De même, les zones où l'on retrouve des courbes fermées à l'intérieur d'une zone chaotique entourant une île sont des zones de résonances (librations). chapitre 1, souschapitre 1, section 5, page 1

Détermination du zéro d'une fonction   (1/2)

Méthodes de la dichotomie et de Newton-Raphson

Elles sont présentées ici. chapitre 1, souschapitre 1, section 5, page 2

Détermination du zéro d'une fonction   (2/2)

Méthode de Newton multidimensionnelle

Méthode de Newton multidimensionnelle

Celle-ci généralise la méthode de Newton-Raphson pour une fonction f : \, R^n \longrightarrow R^p. Le schéma est donné par l'équation suivante : x_{n+1} = x_n - [Df(x_n)]^{-1} \cdot f(x_n), où Df est la matrice jacobienne de f. L'inversion de la matrice jacobienne a été effectuée par décomposition LU. La suite {(x_n)}_n converge vers le zéro de f (sous réserve d'existence et d'unicité).