Résolution de l'équation de Kepler

Restez toujours informé : suivez-nous sur Google (☆)

Introduction

Dans le mouvement keplerien, l'équation de Kepler relie l'anomalie moyenne M = nt à l'anomalie excentrique E par l'équation

M = E - e \cdot sin E

où e est l'excentricité de la planète.

Résoudre cette équation, c'est trouver E(e,M) :

comme série de Fourier puisque c'est une fonction périodique impaire de M
comme série de puissance de e, si e < eo := 0.6627..., rayon de convergence de la série.
comme une valeur numérique avec un nombre de chiffres (d), pour un temps de calcul tc(d) optimisé.

Série de Fourier

C'est Lagrange qui trouve l'expression, bien que le nom J_n(x) soit associé au nom de Bessel.

E-M = fonction impaire périodique de M :

E - M = e \cdot sin E = 2 \cdot Σ_{n = 1} \frac{J _{n} ( n e )}{n} sin (n M)

J_n(x) est la fonction de Bessel de 1ere espèce d'ordre n.

Série entière de l'excentricité

C'est encore Lagrange qui trouve la solution ; et Laplace donnera le rayon de convergence. Ces travaux inspireront Cauchy, qui fondra la théorie des séries analytiques pour résoudre ce problème épineux, celui-ci verra son aboutissement avec les travaux de Puiseux.

L'application du théorème d'inversion de série de Lagrange fournit:

E - M = Σ_{n = 1} \frac{e ^{n}}{n !} \cdot a_{n} (M)

avec $a_{n} (M) = \frac{d ^{n - 1}}{d M ^{n - 1}} (sin^{n} M)$

Le rayon de convergence minimum de la série, qui dépend de M, est atteint pour M= π/2, et vaut: eo = 0.6627434193

tel qu'indiqué par Laplace (1823) et démontré par Cauchy et Puiseux :

et x tel que $\frac{1 - x}{x} = exp (\frac{2}{1 - 2 x})$ .

Ceci rend cette formule inapplicable pour déterminer la position des comètes, dont l'excentricité est souvent voisine de 1.

Les premiers termes sont :

E - M = e \cdot sin M + e^{2} \cdot sin M \cdot cos M + e^{3} \cdot (sin M \cdot cos^{2} M - \frac{sin ^{3} M}{2}) + ...

Note: il est possible d'obtenir ce développement en série en remplaçant dans la série de Fourier précédente les fonctions de Bessel par leur développement limité :

J_{n} (x) = (x /2)^{n} p = 0 \sum \infty \frac{( - 1 ) ^{p}}{p ! ( n + p )!} (x /2)^{2 p}

On obtient alors le développement limité beaucoup plus simplement que par la méthode de l'inversion de série :

E - M = (e - \frac{e ^{3}}{8} + \frac{e ^{5}}{192} + ...) \cdot sin M + (\frac{e ^{2}}{2} - \frac{e ^{4}}{6} + \frac{e ^{6}}{48} + ...) \cdot sin 2 M + (\frac{3}{8} e^{3} - \frac{27}{128} e^{5} + ...) \cdot sin 3 M

+ (\frac{e ^{4}}{3} - \frac{4}{15} e^{6} + ...) \cdot sin 4 M + (\frac{125}{384} e^{5} + ...) \cdot sin 5 M + (\frac{27}{80} e^{6} + ...) \cdot sin 6 M

Il est à noter que bien que la série de Fourier converge pour 0<e<1, et que les développements des fonctions de Bessel aient un rayon de convergence infini, le résultat après réorganisation des termes ne converge que pour e<0.662...

Cas des comètes :
$e > e_{o}$

Le premier à se confronter au problème est Horrocks, puis surtout Halley (1705), pour les calculs sur sa comète d'excentricité e = 0,9673.

Plusieurs solutions ont été proposés en modifiant légèrement la solution de Barker (e = 1). La solution proposée par Bessel (1805) couvre le domaine e > 0.997. Gauss (1809) s'illustra en donnant une belle solution pour 0,2 < e < 0,95

Une généralisation de l'équation de Barker est un développement en série convergeant d'autant plus rapidement que l'excentricité e est proche de 1, ce qui s'avère bien adapté aux cas des comètes (cette série s'applique aussi aux orbites légèrement hyperboliques):

n (t - t_{0}) = S + (1 - 2 \cdot γ) \cdot \frac{S ^{3}}{3} - γ \cdot (2 - 3 \cdot γ) \cdot \frac{S ^{5}}{5} + γ^{2} \cdot (3 - 4 \cdot γ) \cdot \frac{S ^{7}}{7} - ... = S + p = 1 \sum \infty (- γ)^{p - 1} \cdot (p - (p + 1) \cdot γ) \cdot \frac{S ^{2 p + 1}}{2 p + 1}

rayon de convergence: $γ \cdot S^{2} < 1$

avec S = tan(v / 2)

n = \frac{k}{2} \cdot \frac{1 + e}{q ^{3/2}}

γ = \frac{1 - e}{1 + e}

v étant l'anomalie vraie, k la constante gravitationnelle de gauss, e et q étant respectivement l'excentricité et le périastre de l'orbite, t le temps et t₀ étant l'instant du passage au périastre.

Lorsque e =1, la série se réduit à l'équation de Barker

Calcul numérique

L'équation de Kepler peut être résolue à l'aide d'un Algorithme de recherche d'un zéro d'une fonction. Les méthodes du type encadrement, Méthode de dichotomie , Méthode de la fausse position nécessitent un encadrement de départ dans lequel la racine est présente. Du fait de la périodicité et la parité de l'équation de Kepler, il est toujours possible de ramener l'intervalle de départ à [0; π]. Ceci fournit un encadrement de départ pour ces méthode, mais il est aisé d'en trouver de plus fins.

Les méthodes du type point fixe nécessitent une estimation de départ de la racine, le germe de la méthode E₀, pour lancer les calculs : il en existe de multiples dans la littérature, dont la plus simple est E₀=M.

La méthode de type point fixe la plus simple, celle utilisée par Kepler, est :

E_{0} = M

E_{n + 1} = E_{n} + e \cdot sin E_{n} pour n = 0, 1 \dots

converge lentement lorsque e est proche de 1. Il est alors avantageux de lui adjoindre un algorithme d'accélération de la convergence :le Delta-2 de Aitken, par exemple, ou la variante de Steffensen.

L'équation de Kepler se prête particulièrement bien aux algorithmes nécessitant le calcul des dérivées successives élevés, du fait du faible coût en calcul machine nécessaire. En effet :

f (E) = E - e \cdot sin E - M

f^{'} (E) = 1 - e \cdot cos E

f^{''} (E) = e \cdot sin E

f^{'''} (E) = e \cdot cos E

\dots

Les dérivées suivantes se déduisant cycliquement des précédentes. Les variantes de la Méthode de Newton et de Halley d'ordre plus élevés sont donc très performantes dans ce cas. Il est à noter que ces méthodes peuvent dans certains cas avoir des difficultés à converger (e proche de 1 et M proche de 0). Il est préférable dans ces zones soit de proposer une valeur de départ moins grossière (germe de Mikkola), soit de brider les méthodes itératives pour les forcer à converger (modification de Hamming de la méthode de Newton), ou d'utiliser des méthodes itératives à convergence moins locale (Méthode de Laguerre).

Dans le cas des comètes, la résolution de la généralisation quasi-parabolique de l'équation de Barker pose deux problèmes :

le calcul approché de la série qui peut nécessiter un grand nombre de termes, voire être impossible si elle diverge. Il se trouve que cette série se prête particulièrement bien à l'utilisation d'algorithmes d'accélération de la convergence, notamment l'algorithme ε de Peter Wynn, qui non seulement accélère la convergence, mais étend son domaine de convergence. En pratique, les difficultés interviennent lorsque la comète est très loin de son périastre (elle est alors depuis longtemps invisible), ou son excentricité diffère notablement de 1 (dans ce cas,il est plus judicieux de résoudre l'équation de Kepler elliptique ou hyperbolique).
La résolution de l'équation proprement dite. Celle-ci peut être effectuée par les méthodes de type Newton avec :

F (S) = - n (t - t_{0}) + S + p = 1 \sum \infty (- γ)^{p - 1} \cdot (p - (p + 1) \cdot γ) \cdot \frac{S ^{2 p + 1}}{2 p + 1}

en remarquant que la dérivée s'exprime simplement:

$F^{'} (S) = \frac{1 + S ^{2}}{( 1 + γ \cdot S ^{2} ) ^{2}}$ , les dérivées suivantes s'en déduisent facilement.

On pourra choisir comme valeur initiale de l'itération S₀, la solution de l'équation cubique obtenue en retenant les premiers termes (différant légèrement de l'équation de Barker), à l'aide de la Méthode de Cardan

n (t - t_{0}) = S_{0} + (1 - 2 \cdot γ) \cdot \frac{S _{0} ^{3}}{3}

Recherches actuelles

Les calculs via les intégrateur symplectiques exigent de rester toujours en butée du nombre de digits, dans le moindre coût de calcul.
Depuis 300 ans, on cherche la « meilleure » méthode. Elle reste à trouver !

Bien sûr, cela dépend beaucoup du doublet (M,e), M compris entre 0 et π et de e, surtout quand e est voisin de 1.

Nijenhuis (1991) adopte la méthode de Mikkola (1987) qui est la méthode de Newton d'ordre 4, en choisissant « adéquatement » le germe Eo en fonction du doublet (M,e).

Il est clair que dans les calculs numériques, le volume de calculs est essentiel, autant que le nombre de décimales, vu l'instabilité du système solaire évaluée à un coefficient de Liapunov de 10^(t/5Myr). On se heurte à une muraille exponentielle : difficile d'aller plus loin que 25 Myr, même avec un traitement 128 bits.

Ce sont ces calculs (astronomiques... mais informatisés) qui tournent sur les machines de l'IMCCE-Paris. Le calcul de l'ensoleillement terrestre à la latitude 65°Nord, I(65,t) est calculé et on essaie d'en déduire la corrélation avec le climat passé : l'échelle géologique jusqu'au Néogène (25M ans) en est déduite (échelle géologique Gradstein 2004). Prochaine étape prévue : les 65 M ans.

Histoire des sciences

Avant Kepler, l'équation était déjà étudiée pour d'autres motifs :

c'est le problème de la réduction des coordonnées locales aux cordonnées géocentriques : il faut réduire la correction de parallaxe. Habash al Hasib s'y est déjà attaqué.

Avant 1700, il y a déjà beaucoup de tentatives : Kepler naturellement, Curtz (1626), Niele, Boulliau (1645, 1657), Seth Ward (1653), Paganus (1657), Horrebow (1717), Cassini (1669), Newton (1665?), Wren (1658), Wallis (1659), Jeremiah Horrocks (1638)...