Étude dynamique des systèmes multiples de petits
corps

Ajustement par moindres carrés Sommaire

Théorie des moindres carrés

Dans ce chapitre, nous allons nous intéresser `a la fa¸con d’utiliser les moindres carrés qui nous permettra d’ajuster notre mod`ele aux observations. Les coefficients nécessaires sont obtenus en intégrant des équations différentielles appelées équations variationnelles. Nous présenterons ensuite les modifications effectuées pour ajuster les éléments elliptiques initiaux de nos corps `a des observations en coordonnées équatoriales de ces mêmes corps. 4.1 Théorie des moindres carrés Nous allons présenter ici le principe sur lequel repose les moindres carrés et les équations variationnelles.

Méthode Préliminaire

Considérons un mod`ele dynamique `a N-corps, dépendant d’un jeu de constantes c = (−→r0i , −→˙ r0i , pk)i=1,N;k=1,N′. Les param`etres (−→r0i , −→˙ r0i)i=1,N désignent les positions et vitesses initiales des N corps du syst`eme. Les param`etres (pk)k=1,N′ désignent les autres constantes dont peuvent dépendre le mod`ele, comme les masses, les coefficients d’aplatissement, les 4. Ajustement par moindres carres´ 52 directions des pˆoles de rotation . . . On dispose d’un nombre L d’observations du syst`eme, que nous allons noter (−→ r o li)l=1,L;i=1,N ). De même, (−→ r c li)l=1,L;i=1,N ) désignera l’ensemble des positions données par le mod`ele pour ces mêmes instants. Ces L observations nous permettent d’évaluer 3L différences entre notre mod`ele et les observations. Ces différences, que l’on désigne aussi sous le nom ’O-C’ (positions observéespositions calculées), seront notées ∆−→rli et sont égales par définition `a : ∆−→rli = −→ r o li − −→ r c li. Il nous faut alors faire l’hypoth`ese tr`es importante que notre mod`ele décrit parfaitement le syst`eme. Les positions calculées dépendent des conditions initiales et param`etres estimés `a tout instant tl par l’intermédiaire d’une fonction f i telle que −→ r c li = f i (tl , c) (4.1) En raison de l’hypoth`ese faite juste avant, on consid`ere que les positions observées dépendent des conditions initiales réelles et de la valeur réelle des param`etres par l’intermédiaire de la même fonction f. En notant c ′ = (−→r0i ′ , −→˙ r0i ′ , p′ k )i=1,N;k=1,N′ les conditions initiales et param`etres réels, on obtient : −→ r o li = f i (tl , c ′ ) (4.2) Si l’on est proche de la solution réelle du syst`eme, c’est `a dire que les constantes estimées sont proches de leur valeur réelle, on peut se contenter d’un développement de Taylor `a l’ordre 1 de −→ r c li au voisinage de cette solution réelle : −→ r o li = f i (tl , c) + 6N X +N′ k=1 ∂fi ∂c′ k     tl ,ck (c ′ k − ck) + O(c ′ k − ck) 2 (4.3) On obtient alors une expression liant les O-C du syst`emes aux écarts entre les valeurs estimées des constantes et leurs valeurs réelles que nous allons noter ∆ck = ck − c ′ k, que l’on appelle équation de condition : ∆ −→rli ≈ 6N X +N′ k=1 ∂fi ∂ck     tl ,ck (c ′ k − ck) = 6N X +N′ k=1 ∂fi ∂ck     tl ,ck ∆ck (4.4) Ceci correspond `a un syst`eme de 3N ∗ L équations `a 6N + N ′ inconnues. Le but de la méthode des moindres carrés est d’estimer la valeur des param`etres qui permet de minimiser les O-C du syst`eme. Il est évident que dans le cas général, la dépendance entre ∆ck et les O-C n’est pas linéaire. L’une des conditions essentielles `a la validité de cette approximation est de se trouver dans un voisinage de la solution réelle. Notation matricielle Nous pouvons reformuler l’équation 4.4 sous forme matricielle, pour cela, nous allons utiliser les notations suivantes : Y i =   ∆ −→r1i . . . ∆ −→rLi   , Xi =   ∂fi ∂c1     t1,c1 . . . ∂fi ∂c6N+N′     t1,c6N+N′ . . . . . . . . . ∂fi ∂c1     tL,c1 . . . ∂fi ∂c6N+N′     tL,c6N+N′   , A =   ∆c1 . . . ∆c6N+N′   (4.5) 4. Ajustement par moindres carres´ 53 Pour chaque corps i du syst`eme, nous allons alors avoir une équation matricielle du type Y i = X iA. Nous cherchons alors les composantes de la matrice A qui permettent de minimiser celle de la matrice globale Y regroupant toutes les matrices Y i . Résolution Dans le cas idéal : – le développement présenté dans l’équation 4.3 est exact `a l’ordre 1, – la fonction f i représente exactement le comportement du syst`eme, – les résidus Y ne sont pas entˆachés d’erreurs de mesures, il est alors possible de trouver des composantes de A qui rendent nuls les résidus Y . En raison du fait que les hypoth`eses présentées ci-dessus ne sont pas vérifiées, le syst`eme d’équations Y = XA n’admet pas de solution exacte. La résolution approchée de Y = XA nous donnera donc une premi`ere solution approchée. Il faudra donc itérer le processus jusqu’`a la convergence vers une valeur la plus petite possible des résidus. En particulier, les observations sont des mesures dont l’on ne connaˆıt pas a priori la loi des erreurs. Néanmoins, comme les observations sont a priori indépendantes entre elles, nous allons supposer que leur loi des erreurs est gaussienne. Pour trouver la valeur de A qui permet de minimiser les résidus, considérons l’application F : A 7→ kY − XAk 2 (4.6) On cherche un minimum de cette application F, donc une matrice A telle que dF = 0. C’est de ce principe que vient le nom moindres carrés, puisque l’on cherche `a minimiser l’écart quadratique entre les résidus et la matrice XA. La méthode des moindres carrés est l’estimateur le plus efficace dans le cas d’une loi des erreurs gaussienne. Si les colonnes de la matrice X sont linéairement indépendantes, alors le syst`eme dF = 0 admet une unique solution. Par définition : kY − AXk 2 = X L l=1 ∆ −→rli − 6N X +N′ k=1 ∂fi ∂ck     tl ,ck ∆ck !2 = X L l=1 Yl − 6N X +N′ k=1 XlkAl !2 (4.7) Il nous faut alors calculer le gradient de cette expression par rapport aux variables Al ∂ ∂Am kY − AXk 2 = X L l=1 ∂ ∂Am Yl − 6N X +N′ k=1 XlkAl !2 = 2X L l=1 Yl − 6N X +N′ k=1 XlkAl ! (−Xlm) = 2X L l=1 −XlmYl + 6N X +N′ k=1 XlmXlkAl ! (4.8) 4. Ajustement par moindres carres´ 54 Un gradient nul peut donc se réécrire comme un syst`eme de 6N + N ′ équations de la forme : X L l=1 XlmYl − 6N X +N′ k=1 XlmXlkAl ! = 0 (4.9) Ou encore sous la forme matricielle : tXY = tXXA (4.10) La matrice A permettant de minimiser les résidus sera donc la matrice A = ( tXX) −1 tXY . La principale difficulté viendra de l’inversion de la matrice tXX qui peut se révéler impossible, ou du moins numériquement difficile lorsque des observations sont non-indépendantes. Par exemple, des observations obtenues au cours d’une même nuit, pour un objet lointain dont le mouvement est faiblement ou non décelable sur une nuit pour une précision donnée de l’observation, correspondront `a des équations de condition presque identiques. Dans ce cas, même si notre programme parvient `a inverser la matrice tXX, les équations de condition n’étant pas linéairement indépendantes, la correction des param`etres n’aura pas de signification réelle. Une autre source possible d’erreur sur la correction obtenue par les moindres carrés est la corrélation trop forte entre deux param`etres. Dans ce cas, il est impossible de séparer l’influence des deux param`etres et les corrections proposées peuvent ne pas avoir de signification physique. Par exemple, si l’on consid`ere le cas d’un satellite en orbite et doté d’une inclinaison nulle autour d’une plan`ete aplatie. Dans ce cas, il est impossible de séparer l’influence de l’aplatissement dynamique et la masse du corps central. Un autre probl`eme qui risque de survenir est la commensurabilité des différents param`etres `a ajuster. En effet, dans le cas o`u les param`etres `a ajuster ont des valeurs trop dissemblables (par exemple, le demi-grand axe de Pluton a une valeur d’environ 39 UA, alors que la masse de Nix est au maximum de 10−12 masse solaire), l’inversion de A peut alors poser probl`eme en raison de l’incapacité de l’ordinateur `a comparer des nombres si différents en double précision. 

Equations variationnelles

Dans cette partie, nous allons nous intéresser plus particuli`erement `a la mani`ere dont on peut obtenir la valeur des coefficients de la matrice X. En effet, il est nécessaire de pouvoir obtenir, pour tout instant tl et pour tout corps i du syst`eme la valeur des coefficients ∂f i ∂ck     tl ,ck . Pour cela, rappelons que nous avons : d 2−→rli dt2 (tl) = −→F mi (tl , −→rl1, −→˙ rl1, . . . , −→rlN , −→˙ rlN , p) (4.11) La dérivée de l’équation précédente par rapport `a une constante quelconque ck donnera : ∂ ∂ck d 2−→rli dt2 (tl) = ∂ ∂ck −→F mi (tl , −→rl1, −→˙ rl1, . . . , −→rlN , −→˙ rlN , p) = 1 mi « X N j=1 ∂ −→F ∂ −→rlj ∂ −→rlj ∂ck + ∂ −→F ∂ −→˙ rlj ∂ −→˙ rlj ∂ck ! + ∂ −→F ∂p ∂p ∂ck # (4.12) 4. Ajustement par moindres carres´ 55 La constante ck étant indépendante du temps, il est possible d’inverser la dérivation par rapport `a ck et par rapport au temps, on aboutit alors `a d 2 dt2  ∂ −→rli ∂ck  = 1 mi « X N j=1 ∂ −→F ∂ −→rlj ∂ −→rlj ∂ck + ∂ −→F ∂ −→˙ rlj ∂ −→˙ rlj ∂ck ! + ∂ −→F ∂p ∂p ∂ck # (4.13) Qui peut aussi s’écrire sous la forme : d 2 dt2 ∂f i ∂ck     tl ,ck ! = 1 mi « X N j=1 ∂ −→F ∂ −→rlj ∂f j (tl) ∂ck + ∂ −→F ∂ −→˙ rlj ∂ −→˙ rlj ∂ck ! + ∂ −→F ∂p ∂p ∂ck # (4.14) Nous nous retrouvons alors `a intégrer 6N + N ′ syst`emes différentiels d’ordre deux (un pour chaque constante du mod`ele), de 3N équations chacun. Ces équations seront intégrées numériquement en même temps que les équations du mouvement. Les conditions initiales de ces syst`emes seront sous la forme (0, . . ., 0, 1, 0, . . ., 0), o`u le 1 sera au rang k. Dans le cas de notre syst`eme, les forces modélisées ne tiennent pas compte de la vitesse. Par conséquent, l’avant dernier terme de nos équations variationnelles est nul. 

Calcul des équations variationnelles

Dans cette partie, nous allons présenter les équations variationnelles qui nous ont été nécessaires pendant notre ajustement. Nous allons décomposer ces équations suivant les param`etres qui nous intéressent et suivant l’interaction qu’elles représentent. Par la suite, nous allons devoir modifier ces équations pour qu’elles puissent être utilisées dans l’ajustement de notre mod`ele. En effet, notre mod`ele est exprimé en coordonnées cartésiennes, et les constantes dont il dépend sont les positions-vitesses initiales et les param`etres physiques des corps. Or, nos observations sont en coordonnées équatoriales et nous allons devoir ajuster en éléments elliptiques.

Interactions entre les centres de masses

Dans un premier temps, nous allons nous intéresser `a la partie la plus simple des interactions du syst`eme : les interactions entre les centres de masses et leurs dérivées par rapport aux positions et vitesses initiales et aux masses. Ainsi, dans l’équation 4.14, les termes de la forme ∂f j (tl)/∂ck, ∂ −→˙ rlj/∂ck et ∂p/∂ck sont donnés par l’intégration de ces équations au cours du temps dont les conditions initiales seront de la forme (0, . . ., 0, 1, 0, . . ., 0), o`u le 1 sera au rang k, il ne nous reste donc plus qu’`a expliciter les termes de la forme ∂ −→F /∂−→˙ rlj , ∂ −→F /∂−→˙ rlj et ∂ −→F /∂p. Dérivées par rapport aux vitesses Les forces prises en compte ne dépendent pas des vitesses des corps, par conséquent, tous les termes de la forme ∂ −→F /∂−→˙ rlj seront nuls. 4. Ajustement par moindres carres´ 56 Dérivées par rapport aux positions Explicitons maintenant les termes de la forme ∂ −→F /∂−→rlj . Considérons les interactions subies par le corps i sous l’influence des autres objets du syst`eme : F~ i mi = X N j=1 −GMj ~rj − ~ri r 3 ij + X 4 l=1,i6=l −Gml ~rl − ~ri r 3 il (4.15) Une fois projeté sur un axe γ, on obtient : Fiγ mi = X N j=1 −GMj γj − γi r 3 ij + X 4 l=1,i6=l −Gml γl − γi r 3 il (4.16) Nous allons alors nous trouver face `a quatre cas de figures lorsque l’on dérive : 1. on dérive Fiγ par rapport `a la position initiale du même corps i, située sur le même axe γ 2. on dérive Fiγ par rapport `a la position initiale d’un autre corps j, située sur le même axe γ 3. on dérive Fiγ par rapport `a la position initiale du même corps i, située sur un autre axe ζ 4. on dérive Fiγ par rapport `a la position initiale d’un autre corps j, située sur un autre axe ζ Dans le premier cas, nous aurons alors : 1 mi ∂Fiγ ∂γi = X N k=1 GMk r 3 ik  1 − 3(γk − γi) 2 r 2 ik  + X 4 l=1,i6=l Gml r 3 il  1 − 3(γl − γi) 2 r 2 il  (4.17) Dans le deuxi`eme cas : 1 mi ∂Fiγ ∂γj = Gmj r 3 ij  1 − 3(γj − γi) 2 r 2 ij  (4.18) Dans le troisi`eme cas : 1 mi ∂Fiγ ∂ζi = X N k=1 GMk r 3 ik 3(γk − γi)(ζk − ζi) r 2 ik + X 4 l=1,i6=l Gml r 3 il 3(γl − γi)(ζl − ζi) r 2 il (4.19) Dans le quatri`eme cas : 1 mi ∂Fiγ ∂ζj = − Gmj r 3 ij 3(γj − γi)(ζj − ζi) r 2 ij (4.20) 4. Ajustement par moindres carres´ 57 Dérivées par rapport aux masses En ce qui concerne les dérivées par rapport aux masses du syst`eme, la dérivée de F~ i est bien sˆur indépendante de mi . Nous aurons donc, pour la masse mj : 1 mi ∂F~ i ∂mj      γ = −G γj − γi r 3 ij (4.21) 4.2.2 Dérivées des forces d’aplatissement Nous ne prendrons ici en compte que l’effet de l’aplatissement polaire J2 pour des raisons expliquées dans le chapitre suivant. L’expression de cette force ne dépendant pas de la vitesse, nous ne considérerons les dérivées que par rapport aux masses, aux coefficients d’aplatissement, aux positions des corps, et aux coordonnées des pˆoles de rotation, les termes de la forme ∂ −→F /∂−→˙ rlj étant, comme pour les interactions entre les centres de masse, nuls.