|
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é
|
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

dans le champ de gravité de deux masses

et

,

, 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
et
.
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
finalement obtenu s'écrit :
avec comme équations associées :
Les variables
et
se rapportent à la position de la particule,
et
sont leurs variables dites "conjuguées" et sont liées à la vitesse de la particule.
Le Hamiltonien
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 :
où
chapitre 1, souschapitre 1, section 1, page 3
Présentation du problème
(3/3)
Intégrale de Jacobi
L'intégrale de Jacobi
est définie par
, où
est le Hamiltonien du système. C'est donc une intégrale première du mouvement. Cela se traduit dans l'espace des phases
par le fait que le système doive rester sur la variété (de dimension 3)
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

et

, on obtient alors :
Puisque

et

, et que la courbe de vitesse nulle s'écrit

,
nous obtenons :
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 :
Cette expression devant être toujours positive, on voit que pour une valeur fixée de J, le plan (Q
1,Q
2) sera partagé en plusieurs régions:
- Celles où l'équation précédente est négative sont des zones où le mouvement n'est pas possible, ce sont des zones interdites.
- Les points (Q1,Q2) pour lesquels l'expression précédente s'annule forment les séparatrices appelées "courbes de vitesse nulle". Notons que ces séparatrices ne sont pas interdites ! En effet, le fait que la vitesse s'y annule ne signifie pas que l'accélération soit elle aussi nulle. Les trajectoires de la particule peuvent donc passées par un point de ces courbes, où elles sont tangentes.
Bien entendu, ces courbes ne sont visible que dans le repère tournant.
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

. 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:
avec
Les points de Lagrange

et

sont les premiers points où les courbes de vitesse nulle apparaissent. En se plaçant en

, on peut donc déterminer le seuil d'apparition des courbes.

étant de coordonnées

et

, on a

et

égaux à 1 et :
Les courbes de vitesse nulle existent donc pour:
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

,

et

. On utilise bien sûr les valeurs numériques de

calculées lors de la détermination des points fixes,

étant nul pour ces points. Ainsi :
- cas 1 : si
, il n'y a pas de courbes de vitesse nulle à déterminer.
- cas 2 : si
,
,
et
, on se trouve dans le cas où deux zones interdites symétrique se développent autour des points de Lagrange
et
. Il suffit donc de déterminer celles autour de
et la seconde est déduite directement par symétrie.
- cas 3 : si
,
,
et
, ces deux zones se rejoignent autour de
et la courbe de vitesse nulle forme une seule courbe en C fermée.
- cas 4 : si
,
,
et
, la courbe en C se referme en
. Deux courbes de vitesse nulles sont à déterminer.
- cas 5 : enfin si
,
,
et
, la courbe finie par se refermée aussi en
et alors trois courbes sont à déterminer, chacune contenant respectivement
, la première et la seconde masse.
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:
- cas 2 : pour
constant égale à
on effectue la recherche de zéro pour f entre
et
- cas 3 : idem
- cas 4 : pour la première courbe,
constant égale à
, on effectue la recherche de zéro pour f entre
et
; pour la seconde,
constant égale à
, on effectue la recherche de zéro pour f entre
et
- cas 5 : pour la première courbe,
constant égale à
, on effectue la recherche de zéro pour f entre
et
; pour la seconde,
constant égale à
, on effectue la recherche de zéro pour f entre
et
; pour la troisième,
constant égale à
, on effectue la recherche de zéro pour f entre
et
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:
or

et

sont reliés via la différentielle de f :
On déduit de ces deux équations le système suivant:
On a ainsi un système de deux équations différentielles à deux inconnues donnant l'évolution de

et

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 :
- Notez que le système a été écrit de façon à évitez une division par zéro,
et
ne pouvant être nulles simultanément.
- On rappelle que :
d'où
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

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 à

et y celui de

. 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

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 :

.
Un intégrateur numérique est dit symplectique s'il conserve la structure hamiltonienne des
équations, c'est-à-dire les expressions des

et des

.
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

), et sont notés

et

.
On passe alors du point

au point

via la relation :
où les

sont donnés par les expressions suivantes :

,

,
et où les

valent (pour la méthode choisie; réf : "Numerical Hamiltonian Problems", J.M. Sanz-Serna et M.P. Calvo, p. 102) :

,

,
Les

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

ne fait apparaître ni

ni

, et celle définissant

ne fait pas intervenir

. On peut donc d'abord déterminer

, puis

, et enfin

.
Implémentation
Les coefficients

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

est alors la solution (sous réserve d'existence et d'unicité) de

.
On détermine ensuite

comme zéro de la fonction
puis

comme zéro de
On peut alors calculer

.
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é
et pour ordonnée la coordonnée généralisé
(ou, quand le problème permet une formulation hamiltonienne en coordonnées action-angle : (angle-action)),
(j
i) é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
(j
i) é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

. Le schéma est donné par l'équation suivante :
![x_{n+1} = x_n - [Df(x_n)]^{-1} \cdot f(x_n)](eq_tex_trois-corps-cours-2007/equation38.png)
, où

est la matrice jacobienne de

. L'inversion de la matrice jacobienne a été effectuée par
décomposition LU. La suite

converge vers le zéro de

(sous réserve d'existence et d'unicité).