Mémoires extraits des recueils de l’Académie royale de Berlin/Théorie des variations séculaires des éléments des Planètes/1

Joseph-Louis Lagrange

Théorie des variations séculaires des éléments des Planètes. (Première Partie.)

Mémoires extraits des recueils de l’Académie royale de Berlin, Gauthier-Villars, 1868, Œuvres de Lagrange. Tome V (p. 125-207).

THÉORIE
DES VARIATIONS SÉCULAIRES
DES ÉLÉMENTS DES PLANÈTES.

---

PREMIÈRE PARTIE

CONTENANT LES PRINCIPES ET LES FORMULES GÉNÉRALES POUR DÉTERMINER CES VARIATIONS.

---

(Nouveaux Mémoires de l’Académie royale des Sciences et Belles-Lettres
de Berlin, année 1781.)

Si les Planètes étaient simplement attirées par le Soleil, et n’agissaient point les unes sur les autres, elles décriraient autour de cet astre des ellipses invariables suivant les lois de Képler, comme Newton l’a démontré le premier, et une foule d’Auteurs après lui. Mais les observations ont prouvé que le mouvement elliptique des Planètes est sujet à de petites variations, et le calcul a démontré que leur attraction mutuelle peut en être la cause. Ces variations sont de deux espèces les unes périodiques et qui ne dépendent que de la configuration des Planètes entre elles ; celles-ci sont les plus sensibles, et le calcul en a déjà été donné par différents Auteurs ; les autres séculaires et qui paraissent aller toujours en augmentant, ce sont les plus difficiles à déterminer tant par les observations que par la Théorie. Les premières ne dérangent point l’orbite primitive de la Planète ; ce ne sont, pour ainsi dire, que des écarts passagers qu’elle fait dans sa course régulière, et il suffit d’appliquer ces variations au lieu de la Planète calculé par les Tables ordinaires du mouvement elliptique. Il n’en est pas de même des variations séculaires. Ces dernières altèrent les éléments mêmes de l’orbite, c’est-à-dire la position et les dimensions de l’ellipse décrite par la Planète ; et quoique leur effet soit insensible dans un court espace de temps, il peut néanmoins devenir à la longue très-considérable.

C’est une Théorie complète de ces sortes de variations que j’entreprends de donner ; objet qui intéresse également les Astronomes par son utilité pour la perfection des Tables, et les Géomètres par les recherches nouvelles d’analyse auxquelles il donne lieu. Quoique j’aie déjà rempli une partie de cet objet dans les Mémoires sur les équations séculaires des nœuds et des inclinaisons des Planètes^[1], et sur l’altération de leurs mouvements moyens^[2] ; et que j’aie même donné, il y a longtemps, dans les Recherches sur les Satellites de Jupiter^[3], et dans celles sur Jupiter et Saturne^[4], des méthodes et des formules générales pour déterminer ce genre d’inégalités, dont on n’avait encore qu’une connaissance imparfaite et peu exacte ; je crois cependant devoir reprendre cette matière en entier, pour la traiter à fond et d’une manière plus directe et plus rigoureuse que je ne l’ai fait. C’est à quoi sont destinées ces nouvelles Recherches, que je divise en deux Parties. Je donnerai dans la première les formules nécessaires pour déterminer les variations des éléments d’une Planète, réduites à la forme la plus générale et la plus simple ; et j’en ferai, dans la seconde, l’application aux variations séculaires des excentricités, des aphélies, des nœuds et des inclinaisons des six Planètes principales.

section première.

méthode générale pour déterminer les variations des éléments
des planètes, causées par leur action mutuelle.

1. On entend par éléments de l’orbite elliptique d’une Planète la moitié du grand axe de l’ellipse, ou la distance moyenne de la Planète au Soleil ; la position de ce grand axe sur le plan de l’orbite, ou le lieu des apsides ; le rapport de la distance des deux foyers au grand axe, ou l’excentricité l’angle que fait avec l’écliptique le plan de l’orbite, ou son inclinaison ; et l’angle que fait avec une ligne fixe, donnée de position sur l’écliptique, l’intersection de ces deux plans, ou la position de la ligne des nœuds. Ces cinq quantités déterminent complètement la grandeur et la position de l’ellipse ; elles sont par conséquent différentes pour les diverses Planètes ; mais elles demeurent les mêmes pour chaque Planète en particulier, du moins tant qu’on fait abstraction des dérangements qu’elle peut éprouver de la part des autres Planètes. Ainsi les quantités dont nous parlons n’entrent point dans les équations différentielles des orbites des Planètes, parce que ces équations sont générales pour toutes les Planètes ; mais elles entrent ensuite comme constantes arbitraires dans les intégrales de ces équations, c’est-à-dire dans les équations algébriques des orbites.

Pour déterminer l’effet des forces perturbatrices d’une Planète sur ses éléments, il n’y aura donc qu’à traiter les équations différentielles de son orbite comme on ferait si ces forces n’existaient pas, et l’on parviendra ainsi à des équations intégrales semblables à celles de l’orbite non troublée, mais dans lesquelles chaque constante arbitraire se trouvera augmentée d’une quantité variable, provenant des forces perturbatrices, et qui exprimera les dérangements causés par ces forces à l’élément de l’orbite représenté par la même constante. De cette manière l’effet total des perturbations sera renfermé dans les variations des éléments ; et pour avoir la partie séculaire de ces variations, il suffira de rejeter tous les termes qui contiendraient des sinus et cosinus, comme ne pouvant donner que des variations périodiques. Tel est, en général, l’esprit de la méthode que je vais développer et appliquer aux Planètes.

2. Considérons d’abord le mouvement d’un corps mû autour d’un centre fixe en vertu d’une force réciproquement proportionnelle au carré de la distance, et dérangé en même temps par des forces perturbatrices données, et très-petites vis-à-vis de la force principale ; ce qui est le cas de toutes les Planètes.

Soient $x,y,z$ les trois coordonnées rectangles qui déterminent la position du corps à chaque instant, et dont l’origine est supposée dans le centre de la force principale ; nommant $g$ la quantité de cette force à la distance $1,$ et $\rho$ la distance du corps au centre, c’est-à-dire le rayon vecteur de l’orbite, en sorte que

$\rho ={\sqrt {x^{2}+y^{2}+z^{2}}},$

on aura ${\frac {g}{\rho ^{2}}}$ pour l’expression générale de cette force, laquelle étant décomposée suivant les trois coordonnées $x,y,z,$ donnera ces trois-ci

${\frac {gx}{\rho ^{3}}},\quad {\frac {gy}{\rho ^{3}}},\quad {\frac {gz}{\rho ^{3}}}.$

Soient de plus toutes les forces perturbatrices réduites à trois, dirigées suivant les mêmes coordonnées, et représentées par $\mathrm {X,Y,Z} .$

Enfin soit $t$ le temps écoulé depuis une époque donnée, et dont les éléments $dt$ soient pris pour constants.

On aura, par les premiers principes de la Dynamique, ces trois équations différentielles du second ordre

${\begin{aligned}{\frac {d^{2}x}{dt^{2}}}+{\frac {gx}{\rho ^{3}}}+\mathrm {X} =&0,\\{\frac {d^{2}y}{dt^{2}}}+{\frac {gy}{\rho ^{3}}}+\mathrm {Y} =&0,\\{\frac {d^{2}z}{dt^{2}}}+{\frac {gz}{\rho ^{3}}}+\mathrm {Z} =&0,\end{aligned}}$

lesquelles serviront à déterminer le mouvement du corps, en vertu des forces ${\frac {g}{\rho ^{2}}},\ \mathrm {X,Y,Z} .$

3. Il faut maintenant intégrer ces équations de manière que, si l’on faisait abstraction des quantités $\mathrm {X,Y,Z} ,$ il en résultât pour l’orbite des équations algébriques.

Pour cela je fais d’abord ces trois combinaisons, qui servent à chasser les termes divisés par $\rho ^{3},$

${\begin{aligned}{\frac {xd^{2}y-yd^{2}x}{dt^{2}}}=&\mathrm {X} y-\mathrm {Y} x,\\{\frac {xd^{2}z-zd^{2}x}{dt^{2}}}=&\mathrm {X} z-\mathrm {Z} x,\\{\frac {yd^{2}z-zd^{2}y}{dt^{2}}}=&\mathrm {Y} z-\mathrm {Z} y.\end{aligned}}$

Il est visible que les premiers membres de ces équations sont intégrables. Je fais donc

(A)

$\left\{{\begin{aligned}d\mathrm {P} =&(\mathrm {Y} z-\mathrm {Z} y)dt,\\d\mathrm {Q} =&(\mathrm {X} z-\mathrm {Z} x)dt,\\d\mathrm {R} =&(\mathrm {X} y-\mathrm {Y} x)dt\,;\end{aligned}}\right.$

substituant ces valeurs et intégrant, j’aurai

(B)

$\left\{{\begin{aligned}{\frac {xdy-ydx}{dt}}=&\mathrm {R} ,\\{\frac {xdz-zdx}{dt}}=&\mathrm {Q} ,\\{\frac {ydz-zdy}{dt}}=&\mathrm {P} .\end{aligned}}\right.$

Les constantes nécessaires pour compléter ces intégrales sont évidemment contenues dans les quantités $\mathrm {P,Q,R} ,$ lesquelles le deviennent elles-mêmes lorsque les forces perturbatrices $\mathrm {X,Y,Z}$ s’évanouissent.

4. En chassant de ces dernières équations les différences $dx,dy,dz,$ ce qui ne demande que de les ajouter ensemble après les avoir multipliées respectivement par $z,-y,x,$ on aura sur-le-champ cette équation finie

(C)

$\mathrm {P} x-\mathrm {Q} y+\mathrm {R} z=0,$

laquelle sous cette forme appartient à un plan passant par l’origine des coordonnées, et dont la position dépend des quantités $\mathrm {P,Q,R} .$ Car il est facile de démontrer que l’inclinaison de ce plan sur celui des $x$ et $y$ aura pour tangente la quantité $\mathrm {\frac {\sqrt {P^{2}+Q^{2}}}{R}} ,$ et que l’intersection des deux plans fera avec l’axe des $x$ un angle dont la tangente sera $\mathrm {\frac {P}{Q}} .$

Le corps se trouvera donc toujours dans ce plan. Or lorsqu’il n’y a point de forces perturbatrices, et que par conséquent les quantités $\mathrm {P,Q,R}$ sont constantes, la position du plan est invariable, et le corps décrit alors une orbite plane. Mais si ces quantités sont variables, la position du plan doit l’être aussi ; cependant elle peut être censée constante pendant que le corps décrit chaque élément de sa trajectoire. Car les équations (B), étant multipliées respectivement par $dz,-dy,dx,$ donnent celle-ci

$\mathrm {P} dx-\mathrm {Q} dy+\mathrm {R} dz=0,$

qu’on voit n’êtreautre chose que la différentielle de l’équation (C) au plan, en y regardant les quantités $\mathrm {P,Q,R}$ comme constantes.

5. Ainsi le plan dont il s’agit passera par chaque élément de l’orbite du corps, et sera celui dont l’intersection avec le plan de l’écliptique se nomme en Astronomie la ligne des nœuds. De sorte que, si l’on nomme $\theta$ la tangente de l’inclinaison de l’orbite sur le plan de projection que nous supposons être celui des $x$ et $y,$ et $\omega$ l’angle de la ligne des nœuds avec une ligne fixe qui est en même temps l’axe des $x,$ on aura

$\theta =\mathrm {\frac {\sqrt {P^{2}+Q^{2}}}{R}} ,\quad \operatorname {tang} \omega =\mathrm {\frac {P}{Q}} ,$

et de là

$\theta \sin \omega =\mathrm {\frac {P}{R}} ,\quad \theta \cos \omega =\mathrm {\frac {Q}{R}} .$

On connaîtra donc les variations des éléments $\theta$ et $\omega$ par le moyen des formules différentielles (A) ; c’est le chemin que nous avons suivi dans le Mémoire sur les nœuds et les inclinaisons des Planètes, imprimé parmi ceux de l’Académie des Sciences de Paris pour l’année 1774.

6. Pour déterminer les variations des autres éléments, il faut trouver de nouvelles intégrales des équations proposées, en sorte qu’on puisse en déduire l’équation algébrique de l’orbite. À cet effet je remarque que la différentiation des quantités ${\frac {x}{\rho }},{\frac {y}{\rho }},{\frac {z}{\rho }},$ dans lesquelles

$\rho ={\sqrt {x^{2}+y^{2}+z^{2}}},$

donne ces formules

${\begin{aligned}d{\frac {x}{\rho }}=&{\frac {y(ydx-xdy)+z(zdx-xdz)}{\rho ^{3}}},\\d{\frac {y}{\rho }}=&{\frac {x(xdy-ydx)+z(zdy-ydz)}{\rho ^{3}}},\\d{\frac {z}{\rho }}=&{\frac {x(xdz-zdx)+y(ydz-zdy)}{\rho ^{3}}}.\end{aligned}}$

Or, par les équations (B), on a déjà

$ydx-xdy=-\mathrm {R} dt,\quad zdx-xdz=-\mathrm {Q} dt,\quad zdy-ydz=-\mathrm {P} dt\,;$

donc on aura

${\begin{aligned}d{\frac {x}{\rho }}=&-{\frac {\mathrm {R} y+\mathrm {Q} z}{\rho ^{3}}}dt,\\d{\frac {y}{\rho }}=&{\frac {\mathrm {R} x-\mathrm {P} z}{\rho ^{3}}}dt,\\d{\frac {z}{\rho }}=&{\frac {\mathrm {Q} x+\mathrm {P} y}{\rho ^{3}}}dt.\end{aligned}}$

Si maintenant, dans ces équations multipliées par $g,$ on substitue à la place des quantités ${\frac {gx}{\rho ^{3}}},{\frac {gy}{\rho ^{3}}},{\frac {gz}{\rho ^{3}}}$ leurs valeurs tirées des équations différentielles données, savoir (2)

$-{\frac {d^{2}x}{dt^{2}}}-\mathrm {X} ,\quad -{\frac {d^{2}y}{dt^{2}}}-\mathrm {Y} ,\quad -{\frac {d^{2}z}{dt^{2}}}-\mathrm {Z} ,$

on aura les transformées suivantes

${\begin{aligned}gd{\frac {x}{\rho }}=&{\frac {\mathrm {R} d^{2}y+\mathrm {Q} d^{2}z}{dt}}+\mathrm {(RY+QZ)} dt,\\gd{\frac {y}{\rho }}=&{\frac {\mathrm {P} d^{2}z-\mathrm {R} d^{2}x}{dt}}+\mathrm {(PZ-RX)} dt,\\gd{\frac {z}{\rho }}=&-{\frac {\mathrm {Q} d^{2}x+\mathrm {P} d^{2}y}{dt}}-\mathrm {(QX+PY)} dt.\end{aligned}}$

Ces équations sont évidemment intégrables lorsque les forces perturbatrices $\mathrm {X,Y,Z}$ sont nulles, auquel cas les quantités $\mathrm {P,Q,R}$ deviennent constantes (3). Je fais donc, en général,

(D)

$\left\{{\begin{aligned}&\mathrm {(RY+QZ)} dt-{\frac {d\mathrm {R} dy+d\mathrm {Q} dz}{dt}}=d\mathrm {N} ,\\&\mathrm {(P\,Z-RX)} dt-{\frac {d\mathrm {P} dz-d\mathrm {R} dx}{dt}}=d\mathrm {M} ,\\-&\mathrm {(QX+PY)} dt+{\frac {d\mathrm {Q} dx+d\mathrm {P} dy}{dt}}=d\mathrm {L} ,\end{aligned}}\right.$

équations qui étant ajoutées respectivement aux précédéntes les rendent intégrables, en sorte qu’on aura les intégrales suivantes

(E)

$\left\{{\begin{aligned}{\frac {gx}{\rho }}=&{\frac {\mathrm {R} dy+\mathrm {Q} dz}{dt}}+\mathrm {N} ,\\{\frac {gy}{\rho }}=&{\frac {\mathrm {P} dz-\mathrm {R} dx}{dt}}+\mathrm {M} ,\\{\frac {gz}{\rho }}=&-{\frac {\mathrm {Q} dx+\mathrm {P} dy}{dt}}+\mathrm {L} .\end{aligned}}\right.$

Les constantes arbitraires que l’intégration introduit sont renfermées dans les quantités $\mathrm {L,M,N} ,$ lesquelles le deviennent elles-mêmes dans le cas où $\mathrm {X,Y,Z}$ sont nulles.

7. Les intégrales qu’on vient de trouver en donnent immédiatement une algébrique, en éliminant les différences $dx,dy,dz\,;$ et pour cela je les ajoute ensemble après les avoir multipliées respectivement par $x,y,z\,;$ ce qui, à cause de

$\rho ={\sqrt {x^{2}+y^{2}+z^{2}}},$

donne celle-ci

$g\rho =\mathrm {R} {\frac {xdy-ydx}{dt}}+\mathrm {Q} {\frac {xdz-zdx}{dt}}+\mathrm {P} {\frac {ydz-zdy}{dt}}+\mathrm {N} x+\mathrm {M} y+\mathrm {L} z\,;$

laquelle, en vertu des équations (B). devient

(F)

$g\rho =\mathrm {N} x+\mathrm {M} y+\mathrm {L} z+\mathrm {P^{2}+Q^{2}+R^{2}} .$

Cette équation considérée sous cette forme est évidemment du second degré, à cause de

$\rho ={\sqrt {x^{2}+y^{2}+z^{2}}}\,;$

de sorte qu’en la combinant avec l’équation linéaire (C), on aura pour la projection de l’orbite une section conique ; et l’orbite elle-même en sera une aussi, tant que les quantités $\mathrm {P,Q,R,L,M,N}$ seront constantes.

8. Pour déterminer dans ce cas la nature et la position de l’orhite, je commence par remarquer que des six constantes dont nous parlons il n’y en a que cinq d’arbitraires. Car, si l’on multiplie les équations intégrales (E) respectivement par $\mathrm {P,-Q,R} ,$ et qu’ensuite on les ajoute ensemble, on aura, en vertu de l’équation déjà trouvée (C), celle-ci

(G)

$0=\mathrm {NP-MQ+LR} ,$

laquelle exprime un rapport général qui doit subsister entre les six quantités $\mathrm {P,Q,R,N,M,L} ,$ soit qu’elles soient constantes ou non.

Cela posé, je fais, pour abréger,

${\begin{aligned}\Pi ^{2}=&\mathrm {P^{2}+Q^{2}+R^{2}} ,\\\lambda ^{2}=&\mathrm {L^{2}+M^{2}+N^{2}} ,\end{aligned}}$

et, prenant trois autres quantités $\mathrm {A,B,C}$ telles que

$\mathrm {A=PM+QN,\quad B=RN-PL,\quad C=-RM-QL} ,$

je suppose

(H)

$\left\{{\begin{aligned}x=&{\frac {\mathrm {P} }{\Pi }}\zeta +{\frac {\mathrm {N} }{\lambda }}\xi +{\frac {\mathrm {C} }{\Pi \lambda }}\psi ,\\y=&-{\frac {\mathrm {Q} }{\Pi }}\zeta +{\frac {\mathrm {M} }{\lambda }}\xi +{\frac {\mathrm {B} }{\Pi \lambda }}\psi ,\\z=&{\frac {\mathrm {R} }{\Pi }}\zeta +{\frac {\mathrm {L} }{\lambda }}\xi +{\frac {\mathrm {A} }{\Pi \lambda }}\psi ,\end{aligned}}\right.$

$\zeta ,\xi ,\psi$ étant de nouvelles variables ; j’aurai d’abord, en vertu de l’équation de condition (G),

${\begin{aligned}\mathrm {P} x-\mathrm {Q} y+\mathrm {R} z=&\Pi \zeta ,\\\mathrm {N} x+\mathrm {M} y+\mathrm {L} z=&\lambda \xi .\end{aligned}}$

De sorte que les deux équations de l’orbite (C) et (F) se réduiront à cette forme très-simple

$\zeta =0,\quad g\rho =\lambda \xi +\Pi ^{2}.$

Or

$\rho ^{2}=x^{2}+y^{2}+z^{2}\,;$

et si l’on substitue les expressions précédentes de $x,y,z,$ qu’on ait égard à l’équation de condition, et qu’on observe que

${\begin{aligned}\mathrm {A^{2}+B^{2}+C^{2}} =&\mathrm {(PM+QN)^{2}+(RN-PL)^{2}+(RM+QL)^{2}} \\=&\mathrm {\left(P^{2}+Q^{2}+R^{2}\right)\left(L^{2}+M^{2}+N^{2}\right)-(PN-QM+RL)^{2}} \\=&\Pi ^{2}\lambda ^{2},\end{aligned}}$

on aura

$\rho ^{2}=\zeta ^{2}+\xi ^{2}+\psi ^{2}\,;$

d’où l’on peut conclure que les quantités $\zeta ,\xi ,\psi$ sont aussi des coordonnées rectangles, qui répondent aux mêmes points que les coordonnées $x,y,z,$ et qui ont la même origine, mais une position différente.

9. L’équation

$\zeta =0$

fait voir que la courbe est toute dans le plan des coordonnées $\xi$ et $\psi \,;$ et l’équation

$g\rho =\lambda \xi +\Pi ^{2},$

dans laquelle

$\rho ={\sqrt {\xi ^{2}+\xi ^{2}}},$

à cause de $\zeta =0,$ montre que cette courbe est une ellipse dans laquelle ${\frac {\Pi ^{2}}{g}}$ est le paramètre du grand axe, ${\frac {\lambda }{g}}$ est l’excentricité, ou la distance des foyers divisée par le grand axe, $\rho$ est le rayon vecteur partant de l’un des foyers, et $\xi$ l’abscisse prise sur le grand axe depuis le même foyer et dirigée vers l’apside supérieure.

Et pour connaître la position de cet axe relativement au plan de projection, il n’y a qu’à supposer, dans les formules (H), les coordonnées $\psi$ nulles, ce qui, à cause de $\zeta =0,$ donne

$x={\frac {\mathrm {N} }{\lambda }}\xi ,\quad y={\frac {\mathrm {M} }{\lambda }}\xi ,\quad z={\frac {\mathrm {L} }{\lambda }}\xi \,;$

or il est visible que, si l’on nomme $\varphi$ l’angle que la projection de cet axe sur le plan des coordonnées $x,y$ fait avec l’axe des $x,$ et $\eta$ l’angle que le même axe fait avec ce plan, on aura

$\operatorname {tang} \varphi ={\frac {y}{x}},\quad \operatorname {tang} \eta ={\frac {z}{\sqrt {x^{2}+y^{2}}}},$

c’est-à-dire

$\operatorname {tang} \varphi =\mathrm {\frac {M}{N}} ,\quad \operatorname {tang} \eta =\mathrm {\frac {L}{\sqrt {M^{2}+N^{2}}}} ,$

d’où, à cause de

$\lambda =\mathrm {\sqrt {L^{2}+M^{2}+N^{2}}} ,$

on tire

$\sin \eta ={\frac {\mathrm {L} }{\lambda }},\quad \cos \eta \sin \varphi ={\frac {\mathrm {M} }{\lambda }},\quad \cos \eta \cos \varphi ={\frac {\mathrm {N} }{\lambda }}.$

Au reste, si l’on substitue ces valeurs de $\mathrm {L,M,N} ,$ ainsi que celles de $\mathrm {P}$ et $\mathrm {Q}$ tirées des formules du n^o 5, dans l’équation de condition (G), on aura

$\theta \cos \eta (\cos \varphi \sin \omega -\sin \varphi \cos \omega )+\sin \eta =0\,;$

d’où résulte la formule

$\operatorname {tang} \eta =\theta \sin(\varphi -\omega ).$

10. Lorsqu’on veut avoir égard à l’effet des forces perturbatrices, les quantités que nous avons supposées constantes ne le sont plus, et l’orbite telle que nous venons de la déterminer variera d’un instant à l’autre, mais elle pourra néanmoins être prise pour invariable pendant que le corps décrit chacun de ses éléments. En effet les équations (E) du n^o  6 étant multipliées respectivement par $dx,dy,dz,$ et ajoutées ensemble, donnent, à cause de $\rho d\rho =xdx+ydy+zdz,$

$gd\rho =\mathrm {N} dx+\mathrm {M} dy+\mathrm {L} dz\,;$

or cette équation est évidemment la différentielle de l’équation (F) ; en y regardant les quantités $\mathrm {L,M,N,P,Q,R}$ comme constantes.

Ainsi, dans le cas de la variabilité de ces quantités, les deux équations (C) et (F) de l’orbite du corps ont la propriété que ces mêmes quantités y peuvent être regardées comme constantes dans la différentiation de ces équations (4). D’où l’on peut conclure que le corps sera mû à chaque instant comme s’il décrivait réellement l’ellipse déterminée par ces équations ; mais cette ellipse variera continuellement de position et de grandeur, et l’on connaîtra les variations de ses éléments au moyen des formules différentielles (A) et (D) des n^os 3 et 6.

11. Les réductions que nous avons faites ci-dessus (8) étant générales, soit que les éléments de l’orbite soient constants ou non, comme on peut s’en convaincre aisément par la nature de nos formules, il s’ensuit que, si dans les formules (H) on substitue pour $\zeta ,\xi ,\psi$ leurs valeurs tirées des équations

$\zeta =0,\quad g\rho =\lambda \xi +\Pi ^{2},\quad \xi ^{2}+\psi ^{2}=\rho ^{2},$

on aura, en général,

${\begin{aligned}x=&{\frac {\mathrm {N} }{\lambda ^{2}}}\left(g\rho -\Pi ^{2}\right)+{\frac {\mathrm {C} }{\lambda ^{2}}}{\sqrt {-\Delta \rho ^{2}+2g\rho -\Pi ^{2}}},\\y=&{\frac {\mathrm {M} }{\lambda ^{2}}}\left(g\rho -\Pi ^{2}\right)+{\frac {\mathrm {B} }{\lambda ^{2}}}{\sqrt {-\Delta \rho ^{2}+2g\rho -\Pi ^{2}}},\\z=&{\frac {\mathrm {L} }{\lambda ^{2}}}\left(g\rho -\Pi ^{2}\right)+{\frac {\mathrm {A} }{\lambda ^{2}}}{\sqrt {-\Delta \rho ^{2}+2g\rho -\Pi ^{2}}},\end{aligned}}$

en faisant, pour abréger,

$\Delta ={\frac {g^{2}-\lambda ^{2}}{\Pi ^{2}}}.$

Ces expressions de $x,y,z$ en $\rho$ sont les résultats des deux équations (C) et (F) combinées avec la formule

$\rho ^{2}=x^{2}+y^{2}+z^{2}$

($\rho$ étant le rayon vecteur) ; or, comme les quantités $\mathrm {L,M,N,P,Q,R}$ demeurent constantes dans la différentiation de ces équations (numéro précédent), il s’ensuit que, pour avoir les différences $dx,dy,dz,$ il suffira de faire varier dans les expressions précédentes la quantité $p,$ en regardant toutes les autres quantités comme constantes, telles qu’elles le seraient dans le cas de l’orbite invariable.

12. Puisqu’on a déjà $x,y,z$ en $\rho ,$ il ne s’agira plus que d’avoir $\rho$ en $t$ pour connaître le lieu du corps dans un instant quelconque. Pour cela on peut se servir d’une quelconque des intégrales trouvées (B) ou (E) : nous choisirons la première des intégrales (B), savoir

${\frac {xdy-ydx}{dt}}=\mathrm {R} ,$

dont le premier membre, en y substituant les valeurs précédentes de $x,y$ et de leurs différentielles, et remarquant que (8)

${\begin{aligned}\mathrm {NB-MC} =&\mathrm {R\left(M^{2}+N^{2}\right)-L(NP-MQ)} \\=&\mathrm {R\left(L^{2}+M^{2}+N^{2}\right)-L(NP-MQ+LR)} \\=&\mathrm {R} \lambda ^{2},\end{aligned}}$

devient

$-{\frac {\mathrm {R} \rho d\rho }{dt{\sqrt {-\Delta \rho ^{2}+2g\rho -\Pi ^{2}}}}}\,;$

de sorte qu’on aura

$dt=-{\frac {\rho d\rho }{\sqrt {-\Delta \rho ^{2}+2g\rho -\Pi ^{2}}}}.$

Cette équation peut aussi se déduire directement des intégrales (B) et (E) sans aucune substitution auxiliaire. Car :

1^o En ajoutant ensemble les carrés des équations (B), on a

${\frac {(xdy-ydx)^{2}}{dt^{2}}}+{\frac {(xdz-zdx)^{2}}{dt^{2}}}+{\frac {(ydz-zdy)^{2}}{dt^{2}}}$

$={\frac {\left(x^{2}+y^{2}+z^{2}\right)(dx^{2}+dy^{2}+dz^{2})}{dt^{2}}}-{\frac {(xdx+ydy+zdz)^{2}}{dt^{2}}}$

$=\mathrm {P^{2}+Q^{2}+R^{2}} =\Pi ^{2},$

c’est-à-dire

${\frac {\rho ^{2}\left(dx^{2}+dy^{2}+dz^{2}\right)}{dt^{2}}}-{\frac {\rho ^{2}d\rho ^{2}}{dt^{2}}}=\Pi ^{2}\,;$

d’où l’on tire

$dx^{2}+dy^{2}+dz^{2}={\frac {\Pi ^{2}}{\rho ^{2}}}dt^{2}+d\rho ^{2}.$

2^o En ajoutant aussi ensemble les carrés des équations (E), mais après y avoir transposé les termes $\mathrm {N,M,L} ,$ on a

$\left({\frac {gx}{\rho }}-\mathrm {N} \right)^{2}+\left({\frac {gy}{\rho }}-\mathrm {M} \right)^{2}+\left({\frac {gz}{\rho }}-\mathrm {L} \right)^{2}=g^{2}-{\frac {2g}{\rho }}(\mathrm {N} x+\mathrm {M} y+\mathrm {L} z)+\lambda ^{2}$

${\begin{aligned}=&{\frac {(\mathrm {R} dy+\mathrm {Q} dz)^{2}}{dt^{2}}}+{\frac {(\mathrm {P} dz-\mathrm {R} dx)^{2}}{dt^{2}}}+{\frac {(\mathrm {Q} dx+\mathrm {P} dy)^{2}}{dt^{2}}}\\=&{\frac {\Pi ^{2}\left(dx^{2}+dy^{2}+dz^{2}\right)}{dt^{2}}}-{\frac {(\mathrm {P} dx-\mathrm {Q} dy+\mathrm {R} dz)^{2}}{dt^{2}}}\,;\end{aligned}}$

mais on a (4)

$\mathrm {P} dx-\mathrm {Q} dy+\mathrm {R} dz=0,$

et (7)

$\mathrm {N} x+\mathrm {M} y+\mathrm {L} z=g\rho -\Pi ^{2}\,;$

donc

$dx^{2}+dy^{2}+dz^{2}=\left({\frac {2g}{\rho }}-\Delta \right)dt^{2}.$

Égalant ces deux valeurs de $dx^{2}+dy^{2}+dz^{2},$ on en tire comme ci-dessus l’équation

$dt=-{\frac {\rho d\rho }{\sqrt {-\Delta \rho ^{2}+2g\rho -\Pi ^{2}}}}.$

13. Cette équation a donc lieu, en général, soit que les éléments de l’orbite soient invariables ou non ; et elle fait voir que les apsides de l’orbite sont rigoureusement dans les points dont les rayons vecteurs sont les racines de l’équation

$-\Delta \rho ^{2}-2g\rho -\Pi ^{2}=0\,;$

de sorte que ${\frac {g}{\Delta }},$ moitié de la somme des deux racines, sera la distance moyenne. C’est aussi ce qui résulte de l’expression même de $\Delta$ (11), en regardant l’orbite comme une ellipse dont ${\frac {\Pi ^{2}}{\rho }}$ est le paramètre et ${\frac {\lambda }{g}}$ l’excentricité (9).

14. Il ne reste donc plus qu’à intégrer l’équation trouvée, pour en déduire la valeur de $\rho$ en $t\,;$ c’est ce qui est facile dans le cas où les quantités $\Delta$ et $\Pi$ sont constantes. Car la quantité sous le signe peut sue mettre sous cette forme

${\frac {g^{2}}{\Delta }}-\Pi ^{2}-\Delta \left(\rho -{\frac {g}{\Delta }}\right)^{2},$

ou bien ${\biggl (}{\biggr .}$à cause de $\left.\Pi ^{2}={\frac {g^{2}-\lambda ^{2}}{\Delta }}\right)$ sous celle-ci

$\Delta \left[{\frac {\lambda ^{2}}{\Delta ^{2}}}-\left(\rho -{\frac {g}{\Delta }}\right)^{2}\right]\,;$

or on sait que

$-{\frac {d\rho }{\sqrt {{\cfrac {\lambda ^{2}}{\Delta ^{2}}}-\left(\rho -{\cfrac {g}{\Delta }}\right)^{2}}}}$

est l’élément de l’angle dont le cosinus est $\rho -{\frac {g}{\Delta }}$ divisé par ${\frac {\lambda }{\Delta }}\,;$ si donc on nomme $\psi ,$ cet angle, on aura

$\rho ={\frac {g+\lambda \cos \psi }{\Delta }},$

et l’équation dont il s’agit deviendra

$dt={\frac {(g+\lambda \cos \psi )d\psi }{\Delta ^{\frac {3}{2}}}},$

laquelle donne en intégrant

$t+T={\frac {g\psi +\lambda \sin \psi }{\Delta ^{\frac {3}{2}}}},$

$T$ étant la constante arbitraire.

C’est la formule ordinaire qui sert de fondement aux solutions du Problème de Képler, et où l’angle $\psi$ est ce qu’on nomme en Astronomie l’anomalie excentrique.

Lorsque les quantités $\Delta$ et $\Pi$ sont variables, on peut aussi employer la même substitution de

$\rho ={\frac {g+\lambda \cos \psi }{\Delta }}\,;$

et l’on trouvera alors

$dt={\frac {g+\lambda \cos \psi }{\Delta ^{\frac {3}{2}}}}\left(d\psi +{\frac {\cos \psi d\lambda }{\lambda \sin \psi }}\right)+{\frac {(g+\lambda \cos \psi )^{2}d\Delta }{\Delta ^{\frac {5}{2}}\lambda \sin \psi }}\,;$

mais cette formule est peu commode pour le calcul, à cause qu’elle contient $\sin \psi$ au dénominateur.

On pourrait remédier à cet inconvénient en substituant directement la valeur de $\rho$ dans les expressions de $x,y$ du n^o 11, ce qui donnera

${\begin{aligned}x=&{\frac {\mathrm {N} }{\Delta \lambda }}(\lambda +g\cos \psi )+{\frac {\mathrm {C} }{{\sqrt {\Delta }}.\lambda }}\sin \psi ,\\y=&{\frac {\mathrm {M} }{\Delta \lambda }}(\lambda +g\cos \psi )+{\frac {\mathrm {B} }{{\sqrt {\Delta }}.\lambda }}\sin \psi ,\end{aligned}}$

et mettant ensuite ces valeurs de $x$ et $y$ dans l’équation

${\frac {xdy-ydx}{dt}}=\mathrm {R} .$

Mais il sera beaucoup plus simple, surtout dans le cas où l’excentricité et l’inclinaison sont fort petites (qui est celui de toutes les Planètes de notre système), d’employer d’abord dans cette même équation la substitution de

$x=r\cos q,\quad y=r\sin q,$

$r$ étant le rayon vecteur de l’orbite projetée et $q$ l’angle décrit par ce rayon ; car on aura par là

$dt={\frac {r^{2}dq}{\mathrm {R} }},$

et il n’y aura plus qu’à mettre pour $r^{2}$ sa valeur en $q$ tirée des deux équations de l’orbite (C) et (F). C’est ainsi que nous en userons dans la suite.

15. Pour appliquer aux Planètes les formules générales que nous venons de trouver, il y faudra substituer les valeurs des forces perturbatrices $\mathrm {X,Y,Z}$ qui résultent de l’attraction que chaque Planète éprouve de la part de toutes les autres.

Soit $\mathrm {S}$ la masse du Soleil, $\mathrm {T}$ celle de la Planète dont on cherche le mouvement, $\mathrm {T',T''} ,\ldots$ les masses des Planètes perturbatrices ; on sait que la Planète $\mathrm {T}$ est attirée vers le Soleil par une force égale à $\mathrm {\frac {S+T}{\rho ^{2}}} ,$ $\rho$ étant sa distance au Soleil, et qu’en vertu de cette force elle doit décrire autour du Soleil la même orbite que si le Soleil était immobile. On peut donc regarder le Soleil comme fixe par rapport à la Planète $\mathrm {T} ,$ mais il faut alors tenir compte de l’action des autres Planètes $\mathrm {T',T''} ,\ldots$ sur le Soleil, en transportant l’effet de cette action à la Planète $\mathrm {T}$ en sens contraire.

Ainsi, prenant le centre du Soleil pour l’origine des coordonnées, et nommant $x,y,z$ celles de l’orbite de la Planète $\mathrm {T}$ autour du Soleil, on aura d’abord dans les formules du n^o 2

$g=\mathrm {S+T} .$

Ensuite, si l’on marque d’un trait toutes les quantités qui se rapportent à la planète $\mathrm {T} ',$ de deux traits celles qui se rapportent à la Planète $\mathrm {T} '',\ldots \,;$ qu’enfin on désigne par $\sigma '$ la distance rectiligne entre les corps $\mathrm {T}$ et $\mathrm {T} ',$ par $\sigma ''$ la distance rectiligne entre les corps $\mathrm {T}$ et $\mathrm {T} '',$ et ainsi du reste ; on trouvera :

1^o Que la force ${\frac {\mathrm {T} '}{\sigma '^{2}}},$ avec laquelle le corps $\mathrm {T} '$ attire le corps $\mathrm {T}$ suivant la direction de la ligne $\sigma ',$ produira ces trois forces suivant les directions des coordonnées $x,y,z,$ savoir

${\frac {\mathrm {T} '(x-x')}{\sigma '^{3}}},\quad {\frac {\mathrm {T} '(y-y')}{\sigma '^{3}}},\quad {\frac {\mathrm {T} '(z-z')}{\sigma '^{3}}}\,;$

2^o Que la force ${\frac {\mathrm {T} '}{\rho '^{2}}},$ avec laquelle la Planète $\mathrm {T} '$ attire le Soleil $\mathrm {S} ,$ étant transportée en sens contraire à la Planète $\mathrm {T} ',$ donnera encore ces trois autres forces suivant les mêmes directions, savoir

${\frac {\mathrm {T} 'x'}{\rho '^{3}}},\quad {\frac {\mathrm {T} 'y'}{\rho '^{3}}},\quad {\frac {\mathrm {T} 'z'}{\rho '^{3}}}.$

On trouvera de pareilles formules pour les forces résultantes de l’attraction des autres Planètes $\mathrm {T'',T'''} ,\ldots \,;$ et, rassemblant respectivement ces différentes forces, on aura les valeurs des forces perturbatrices $\mathrm {X,Y,Z}$ de la Planète $\mathrm {T} ,$ lesquelles seront donc exprimées ainsi

${\begin{aligned}\mathrm {X} =&\mathrm {T} '\left({\frac {x-x'}{\sigma '^{3}}}+{\frac {x'}{\rho '^{3}}}\right)+\mathrm {T} ''\left({\frac {x-x''}{\sigma ''^{3}}}+{\frac {x''}{\rho ''^{3}}}\right)+\mathrm {T} '''\left({\frac {x-x'''}{\sigma '''^{3}}}+{\frac {x'''}{\rho '''^{3}}}\right)+\ldots ,\\\mathrm {Y} =&\mathrm {T} '\left({\frac {y-y'}{\sigma '^{3}}}\,+{\frac {y'}{\rho '^{3}}}\right)\,+\mathrm {T} ''\left({\frac {y-y''}{\sigma ''^{3}}}+{\frac {y''}{\rho ''^{3}}}\right)+\mathrm {T} '''\left({\frac {y-y'''}{\sigma '''^{3}}}+{\frac {y'''}{\rho '''^{3}}}\right)+\ldots ,\\\mathrm {Z} =&\mathrm {T} '\left({\frac {z-z'}{\sigma '^{3}}}\,+{\frac {z'}{\rho '^{3}}}\right)\,+\mathrm {T} ''\left({\frac {z-z''}{\sigma ''^{3}}}\,+{\frac {z''}{\rho ''^{3}}}\right)+\mathrm {T} '''\left({\frac {z-z'''}{\sigma '''^{3}}}+{\frac {z'''}{\rho '''^{3}}}\right)+\ldots ,\\\end{aligned}}$

16. Or on a

${\begin{aligned}\rho \ \ &={\sqrt {x^{2}\,+y^{2}\,+z^{2}}},\\\rho '\ &={\sqrt {x'^{2}+y'^{2}+z'^{2}}},\\\ldots &\ldots \ldots \ldots \ldots \ldots \ldots ,\\\sigma '\ &={\sqrt {(x-x'\,)^{2}+(y-y'\,)^{2}+(z-z'\ )^{2}}},\\\sigma ''&={\sqrt {(x-x'')^{2}+(y-y'')^{2}+(z-z'')^{2}}},\\\ldots &\ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots .\end{aligned}}$

Donc, si l’on fait

$\Omega =\mathrm {T} '\left({\frac {xx'+yy'+zz'}{\rho '^{3}}}-{\frac {1}{\sigma '}}\right)+\mathrm {T} ''\left({\frac {xx''+yy''+zz''}{\rho ''^{3}}}-{\frac {1}{\sigma ''}}\right)$

$+\mathrm {T} '''\left({\frac {xx'''+yy'''+zz'''}{\rho '''^{3}}}-{\frac {1}{\sigma '''}}\right)+\ldots ,$

et qu’on dénote à l’ordinaire par ${\frac {d\Omega }{dx}},{\frac {d\Omega }{dy}},{\frac {d\Omega }{dz}}$ les coefficients de $dx,dy,dz,$ dans la différentielle de la quantité $\Omega$ regardée comme fonction des variables $x,y,z\,;$ il est clair que les expressions précédentes de $\mathrm {X,Y,Z}$ se réduiront à celles-ci

$\mathrm {X} ={\frac {d\Omega }{dx}},\quad \mathrm {Y} ={\frac {d\Omega }{dy}},\quad \mathrm {Z} ={\frac {d\Omega }{dz}}.$

17. On substituera donc ces valeurs dans les équations différentielles (A) et (D), et l’on aura en premier lieu

${\begin{aligned}d\mathrm {P} =&\left({\frac {d\Omega }{dy}}z-{\frac {d\Omega }{dz}}y\right)dt,\\d\mathrm {Q} =&\left({\frac {d\Omega }{dx}}z-{\frac {d\Omega }{dz}}x\right)dt,\\d\mathrm {R} =&\left({\frac {d\Omega }{dx}}y-{\frac {d\Omega }{dy}}x\right)dt.\end{aligned}}$

Ces formules serviront à déterminer les variations de la longitude $\omega$ du nœud et de la tangente $\theta$ de l’inclinaison, en faisant (5)

$\mathrm {P=R} \theta \sin \omega ,\quad \mathrm {Q=R} \theta \cos \omega ,$

et de plus la variation du paramètre ${\frac {\Pi ^{2}}{g}}$ de l’orbite (9), $\Pi ^{2}$ étant égal à $\mathrm {P^{2}+Q^{2}+R^{2}} .$

18. Si l’on voulait déterminer directement les variations de $\Pi ,$ il n’y aurait qu’à ajouter ensemble les équations précédentes, après les, avoir multipliées respectivement par $\mathrm {P,Q,R} .$ On aura ainsi en ordonnant les termes

$\Pi d\Pi ={\frac {d\Omega }{dx}}(\mathrm {Q} z+\mathrm {R} y)dt+{\frac {d\Omega }{dy}}(\mathrm {P} z-\mathrm {R} x)dt-{\frac {d\Omega }{dz}}(\mathrm {P} y+\mathrm {Q} x)dt.$

Or, en substituant pour $\mathrm {P,Q,R}$ les valeurs données par les équa-

tions(B) du n^o 3, on a

${\begin{aligned}(\mathrm {Q} z+\mathrm {R} y)dt=&(xdz-zdx)z+(xdy-ydx)y\\=&x(ydy+zdz)-dx\left(y^{2}+z^{2}\right)=x\rho d\rho -\rho ^{2}dx,\\(\mathrm {P} z-\mathrm {R} x)dt=&(ydz-zdy)z-(xdy-ydx)x\\=&y(xdx+zdz)-dy\left(x^{2}+z^{2}\right)=y\rho d\rho -\rho ^{2}dy,\\(\mathrm {P} y+\mathrm {Q} x)dt=&(ydz-zdy)y+(xdz-zdx)x\\=&-z(xdx+ydy)+dz\left(x^{2}+y^{2}\right)=-z\rho d\rho +\rho ^{2}dz.\end{aligned}}$

Donc si l’on fait, pour abréger,

${\frac {d\Omega }{dx}}x+{\frac {d\Omega }{dy}}y+{\frac {d\Omega }{dz}}z=\Phi ,\quad {\frac {d\Omega }{dx}}dx+{\frac {d\Omega }{dy}}dy+{\frac {d\Omega }{dz}}dz=(d\Omega )$

[l’expression $(d\Omega )$ indique la différentielle partielle de $\Omega ,$ en n’y faisant varier que les $x,y,z$ relatives à la Planète $\mathrm {T}$ ], on aura

$\Pi d\Pi =\Phi \rho d\rho -\rho ^{2}(d\Omega ).$

19. En second lieu, les équations (D) deviendront, par la substitution des valeurs précédentes de $d\mathrm {P} ,d\mathrm {Q} ,d\mathrm {R} ,$

${\begin{alignedat}{4}-{\frac {d\Omega }{dx}}(ydy+zdz)+&{\frac {d\Omega }{dy}}(xdy+\mathrm {R} dt)&+&{\frac {d\Omega }{dz}}(xdz+\mathrm {Q} dt)&=&d\mathrm {N} ,\\{\frac {d\Omega }{dx}}(ydx-\mathrm {R} dt)-&{\frac {d\Omega }{dy}}(xdx+zdz)&+&{\frac {d\Omega }{dz}}(ydz+\mathrm {P} dt)&=&d\mathrm {M} ,\\{\frac {d\Omega }{dx}}(zdx-\mathrm {Q} dt)+&{\frac {d\Omega }{dy}}(zdy-\mathrm {P} dt)&-&{\frac {d\Omega }{dz}}(ydy+xdx)&=&d\mathrm {L} ,\end{alignedat}}$

et si l’on y substitue encore les valeurs de $\mathrm {P} dt,\mathrm {Q} dt,\mathrm {R} dt$ (3), elles se réduiront à la forme suivante

${\begin{aligned}d\mathrm {N} =&2x(d\Omega )-\Phi dx-{\frac {d\Omega }{dx}}\rho d\rho ,\\d\mathrm {M} =&2y(d\Omega )-\Phi dy-{\frac {d\Omega }{dy}}\rho d\rho ,\\d\mathrm {L} =&2z(d\Omega )-\Phi dz-{\frac {d\Omega }{dz}}\rho d\rho ,\end{aligned}}$

les quantités $\Phi$ et $(d\Omega )$ étant les mêmes que dans le numéro précédent.

Par ces équations, on aura les variations de l’excentricité ${\frac {\lambda }{g}},$ de la longitude de l’aphélie $\varphi ,$ et de la latitude $\eta$ de cet aphélie par rapport au plan de projection, au moyen des formules (10)

$\mathrm {L} =\lambda \sin \eta ,\quad \mathrm {M} =\lambda \cos \eta \sin \varphi ,\quad \mathrm {N} =\lambda \cos \eta \cos \varphi .$

20. Au reste il n’est pas nécessaire d’employer ces trois équations pour la détermination des éléments dont il s’agit ; car nous avons vu (8) qu’il y a entre les quantités $\mathrm {L,M,N}$ cette équation de condition

$\mathrm {NP-MQ+LR} =0,$

au moyen de laquelle on peut déterminer une quelconque d’entre elles par les deux autres. On pourrait peut-être douter si en effet les formules différentielles trouvées satisfont à cette équation ; pour lever ce doute et confirmer en même temps la justesse de nos formules, nous allons faire voir comment elles y satisfont ; à cet effet nous ajouterons d’abord ensemble les trois formules du numéro précédent, après les avoir multipliées respectivement par $\mathrm {P,-Q,R} \,;$ ce qui, en ayant égard (4) aux équations

$\mathrm {P} x-\mathrm {Q} y+\mathrm {R} z=0,\quad \mathrm {P} dx-\mathrm {Q} dy+\mathrm {R} dz=0,$

donnera

$\mathrm {P} d\mathrm {N} -\mathrm {Q} d\mathrm {M} +\mathrm {R} d\mathrm {L} =-\left(\mathrm {P} {\frac {d\Omega }{dx}}-\mathrm {Q} {\frac {d\Omega }{dy}}+\mathrm {R} {\frac {d\Omega }{dz}}\right)\rho d\rho \,;$

ensuite les trois formules du n^o 17, étant multipliées respectivement par $\mathrm {N,-M,L}$ et ajoutées ensemble, donneront

$\mathrm {N} d\mathrm {P} -\mathrm {M} d\mathrm {Q} +\mathrm {L} d\mathrm {R} =$

${\frac {d\Omega }{dx}}(\mathrm {L} y-\mathrm {M} z)dt+{\frac {d\Omega }{dy}}(\mathrm {N} z-\mathrm {L} x)dt+{\frac {d\Omega }{dz}}(\mathrm {M} x-\mathrm {N} y)dt\,;$

mais, en substituant les valeurs de $\mathrm {L,M,N}$ résultantes des formules (E) du n^o 6, on a

${\begin{aligned}(\mathrm {L} \,y-\mathrm {M} z)dt=&(\mathrm {Q} dx+\mathrm {P} dy)y-(\mathrm {R} dx+\mathrm {P} dz)z,\\(\mathrm {N} \,z-\mathrm {L} \,x)dt=&-(\mathrm {R} dy+\mathrm {Q} dz)z-(\mathrm {Q} dx+\mathrm {P} dy)x,\\(\mathrm {M} x-\mathrm {N} y)dt=&(\mathrm {R} dx-\mathrm {P} dz)x+(\mathrm {R} dy+\mathrm {Q} dz)y,\end{aligned}}$

ce qui, à cause de

$\mathrm {P} x-\mathrm {Q} y+\mathrm {R} z=0,$

se réduit à

$(\mathrm {L} y-\mathrm {M} z)dt=\mathrm {P} \rho dr,\quad (\mathrm {N} z-\mathrm {L} x)dt=-\mathrm {Q} \rho dr,\quad (\mathrm {M} x-\mathrm {N} y)dt=\mathrm {R} \rho dr\,;$

de sorte qu’on aura

$\mathrm {N} d\mathrm {P} -\mathrm {M} d\mathrm {Q} +\mathrm {L} d\mathrm {R} =\left(\mathrm {P} {\frac {d\Omega }{dx}}-\mathrm {Q} {\frac {d\Omega }{dy}}+\mathrm {R} {\frac {d\Omega }{dz}}\right)\rho d\rho .$

Donc enfin

$\mathrm {N} d\mathrm {P} -\mathrm {M} d\mathrm {Q} +\mathrm {L} d\mathrm {R} +\mathrm {P} d\mathrm {N} -\mathrm {Q} d\mathrm {M} +\mathrm {R} d\mathrm {L} =0,$

équation qui est, comme on voit, la différentielle de celle qu’il s’agissait de vérifier.

21. Puisque

$\lambda ^{2}=\mathrm {L^{2}+M^{2}+N^{2}} ,$

on aura, en ajoutant ensemble les trois équations du n^o 19, après les avoir respectivement multipliées par $\mathrm {N,M,L} ,$

${\begin{aligned}\lambda d\lambda =&2(\mathrm {N} x+\mathrm {M} y+\mathrm {L} z)(d\Omega )-\Phi (\mathrm {N} dx+\mathrm {M} dy+\mathrm {L} dz)\\&-\left(\mathrm {N} {\frac {d\Omega }{dx}}+\mathrm {M} {\frac {d\Omega }{dy}}+\mathrm {L} {\frac {d\Omega }{dz}}\right)\rho dr\,;\end{aligned}}$

or l’équation de l’orbite donne (7, 10)

$\mathrm {N} x+\mathrm {M} y+\mathrm {L} z=g\rho -\Pi ^{2},\quad \mathrm {N} dx+\mathrm {M} dy+\mathrm {L} dz=g\rho d\rho \,;$

de plus, les équations (E) donnent, en substituant pour $\mathrm {P.Q,R}$ leurs valeurs tirées des équations (B),

(I)

$\left\{{\begin{aligned}\mathrm {N} =&x\left({\frac {g}{\rho }}-{\frac {dx^{2}+dy^{2}+dz^{2}}{dt^{2}}}\right)+{\frac {\rho d\rho dx}{dt^{2}}},\\\mathrm {M} =&\,y\left({\frac {g}{\rho }}-{\frac {dx^{2}+dy^{2}+dz^{2}}{dt^{2}}}\right)+{\frac {\rho d\rho dy}{dt^{2}}},\\\mathrm {L} =&\ z\left({\frac {g}{\rho }}-{\frac {dx^{2}+dy^{2}+dz^{2}}{dt^{2}}}\right)+{\frac {\rho d\rho dz}{dt^{2}}},\end{aligned}}\right.$

et les mêmes équations donnent aussi, comme on l’a vu dans le n^o 12,

${\frac {dx^{2}+dy^{2}+dz^{2}}{dt^{2}}}={\frac {2g}{\rho }}-\Delta \,;$

donc, faisant ces substitutions, on aura

$\lambda d\lambda =\left(2g\rho -2\Pi ^{2}-{\frac {\rho ^{2}d\rho ^{2}}{dt^{2}}}\right)(d\Omega )-\Delta \Phi \rho d\rho \,;$

mais on a par les formules du même numéro cité

${\frac {\rho ^{2}d\rho ^{2}}{dt^{2}}}=-\Delta \rho ^{2}+2g\rho -\Pi ^{2}\,;$

donc

$\lambda d\lambda =\left(\Delta \rho ^{2}-\Pi ^{2}\right)(d\Omega )-\Delta \Phi \rho d\rho ,$

$\Delta$ étant égal à ${\frac {g^{2}-\lambda ^{2}}{\Pi ^{2}}}$ (11). Ainsi l’on aura directement par cette équation la variation de l’excentricité ${\frac {\lambda }{g}}.$

22. Enfin, puisque

$\lambda ^{2}=g^{2}-\Delta \Pi ^{2},$

on aura, en substituant la valeur de $\Pi d\Pi$ du n^o 18,

$\lambda d\lambda =-\Delta \Pi d\Pi -\Pi ^{2}{\frac {d\Delta }{2}}=-\Delta \Phi \rho d\rho +\Delta \rho ^{2}(d\Omega )-\Pi ^{2}{\frac {d\Delta }{2}}\,;$

et, cette valeur de $\lambda d\lambda$ étant comparée à celle que nous venons de trouver dans le numéro précédent, il viendra cette formule très-simple

${\frac {d\Delta }{2}}=(d\Omega ),$

laquelle servira à déterminer les variations de la distance moyenne ${\frac {g}{\Delta }}$ (13).

23. La méthode que j’ai suivie pour trouver les formules différentielles des variations des éléments des Planètes est, ce me semble, la plus directe et la plus naturelle qu’il est possible, étant déduite des principes mêmes de la chose ; mais j’aurais pu y parvenir plus simplement par la méthode générale dont je me suis servi pour déterminer les variations de la distance moyenne, dans les Mémoires de 1776^[5], méthode qui a l’avantage d’être applicable à toutes les questions du même genre.

Suivant cette méthode, si $\mathrm {V} =a$ est une intégrale quelconque des équations différentio-différentielles de l’orbite non troublée, $\mathrm {V}$ étant une fonction connue de $t,x,y,z$ et de ${\frac {dx}{dt}},{\frac {dy}{dt}},{\frac {dz}{dt}},$ et $a$ une constante arbitraire, et qu’on veuille supposer $a$ variable pour l’orbite troublée, il n’y aura qu’à différentier, en faisant varier d’un côté la quantité $a,$ et de l’autre les quantités ${\frac {dx}{dt}},{\frac {dy}{dt}},{\frac {dz}{dt}}$ contenues dans $\mathrm {V} ,$ et substituer ensuite à la place des différences de celles-ci les valeurs de ${\frac {d^{2}x}{dt}},{\frac {d^{2}y}{dt}},{\frac {d^{2}z}{dt}},$ dues uniquement aux forces perturbatrices. Car la différentielle de $\mathrm {V} ,$ prise en faisant tout varier, doit devenir nulle d’elle-même, lorsqu’on y substitue pour $d{\frac {dx}{dt}},d{\frac {dy}{dt}},d{\frac {dz}{dt}}$ les valeurs tirées des équations de l’orbite non troublée, puisque $\mathrm {V} =a$ est (hypothèse) une intégrale de ces équations. De sorte qu’en ajoutant à ces valeurs les termes $-\mathrm {X} dt,$$-\mathrm {Y} dy,-\mathrm {Z} dt$ dus aux forces perturbatrices (2), on aura simplement

$d\mathrm {V} ={\frac {d\mathrm {V} }{d{\cfrac {dx}{dt}}}}(-\mathrm {X} dt)+{\frac {d\mathrm {V} }{d{\cfrac {dy}{dt}}}}(-\mathrm {Y} dt)+{\frac {d\mathrm {V} }{d{\cfrac {dz}{dt}}}}(-\mathrm {Z} dt)\,;$

les quantités ${\frac {d\mathrm {V} }{d{\cfrac {dx}{dt}}}},\ldots$ exprimant suivant la notation reçue tes coeffieients de $d{\frac {dx}{dt}},\ldots$ dans la différentielle de $\mathrm {V} .$

De cette manière, si l’on change, suivant le n^o 16, les quantités $\mathrm {X,Y,Z}$ en ${\frac {d\Omega }{dx}},{\frac {d\Omega }{dy}},{\frac {d\Omega }{dz}},$ on aura, pour la variation de $a,$ cette équation

$da=-{\frac {d\mathrm {V} }{d{\cfrac {dx}{dt}}}}{\frac {d\Omega }{dx}}dt-{\frac {d\mathrm {V} }{d{\cfrac {dy}{dt}}}}{\frac {d\Omega }{dy}}dt-{\frac {d\mathrm {V} }{d{\cfrac {dz}{dt}}}}{\frac {d\Omega }{dz}}dt.$

Or on sait que les éléments de l’orbite ne sont autre chose que les constantes arbitraires introduites par l’intégration des équations différentielles primitives, ou des fonctions de ces constantes ; on aura donc les formules différentielles de la variation des éléments, pourvu qu’on ait les valeurs de ces éléments exprimées par les variables de l’orbite non troublée et par leurs différences premières ; valeurs qui peuvent se trouver immédiatement par l’intégration même des équations différentielles de l’orbite, ou bien se déduire des équations finies de l’orbite combinées avec les différences premières de ces équations.

Ayant déterminé ainsi la valeur variable de chaque constante ou élément de l’orbite, on aura, entre ces constantes et $t,x,y,z,{\frac {dx}{dt}},{\frac {dy}{dt}},{\frac {dz}{dt}},$ des équations de la même forme pour l’orbite troublée que pour l’orbite non troublée ; de sorte que les équations finies déduites de eelle-ci, ainsi que les différences premières de ces équations seront encore de la même forme dans les deux cas ; par conséquent les constantes dont nous parlons pourront toujours être regardées et traitées comme invariables dans la différentiation des équations finies de l’orbite ; ce qui rend raison de la remarque faite dans les n^os 4 et 10 sur les équations (C) et (F).

24. Cela posé, si dans les équations différentielles du n^o  3 on suppose $\mathrm {X,Y,Z}$ nulles, on a celles de l’orbite non troublée ; et, intégrant ces équations par la méthode des n^os 3 et 6, on aura des intégrales de la forme (B) et (C), dans lesquelles $\mathrm {P,Q,R,L,M,N}$ seront des constantes arbitraires et invariables. Or les intégrales (B) ont déjà la forme demandée $\mathrm {V} =a\,;$ donc il n’y aura, suivant la méthode précédente, qu’à les différentier, en y faisant varier $\mathrm {P,Q,R}$ et ${\frac {dx}{dt}},{\frac {dy}{dt}},{\frac {dz}{dt}},$ et substituer ensuite $-{\frac {d\Omega }{dx}}dt,\ -{\frac {d\Omega }{dy}}dt,\ -{\frac {d\Omega }{dz}}dt$ à la place de $d{\frac {dx}{dt}},\ d{\frac {dy}{dt}},\ d{\frac {dz}{dt}}\,;$ on aura ainsi sur-le-champ les formules différentielles du n^o  17.

À l’égard des intégrales (E), pour les réduire à la forme dont il s’agit, il ne faut qu’y substituer les valeurs de $\mathrm {P,Q,R}$ données par les intégrales précédentes ; on aura ainsi les équations (I) du n^o  21, lesquelles étant différentiées en y faisant varier $\mathrm {L,M,N}$ et ${\frac {dx}{dt}},\ {\frac {dy}{dt}},\ {\frac {dz}{dt}},$ donneront immédiatement les formules différentielles du n^o 19, en se souvenant que

$\rho dr=xdx+yxy+zdz.$

On peut de même déterminer directement les variations des éléments $\Pi$ et $\Delta$ par le moyen des formules trouvées dans le n^o 12. Car

1^o On a

$\Pi ^{2}={\frac {\rho ^{2}\left(dx^{2}+dy^{2}+dz^{2}\right)}{dt^{2}}}-{\frac {(\rho d\rho )^{2}}{dt^{2}}},$

et cette équation, en y faisant varier $\Pi$ de $d\Pi$ et ${\frac {dx}{dt}},\ {\frac {dy}{dt}},\ {\frac {dz}{dt}}$ comme ci-dessus, donne immédiatement la formule différentielle trouvée à la fin du n^o 18 ;

2^o On a

$\Delta ={\frac {2g}{\rho }}-{\frac {dx^{2}+dy^{2}+dz^{2}}{dt^{2}}},$

et il est visible qu’il en naîtra sur-le-champ la formule différentielle du n^o 22, par la variation de $\Delta ,\ {\frac {dx}{dt}},\ {\frac {dy}{dt}},\ {\frac {dz}{dt}}.$

25. Il n’est pas nécessaire au reste, pour l’usage de la méthode précédente, que les intégrales soient réduites à la forme $\mathrm {V} =a,$ ainsi que nous l’avons supposé ; mais cette réduction, qui est d’ailleurs toujours possiblé, sert à rendre le calcul plus direct et les formules plus simples. Comme cette méthode peut être d’une grande utilité dans plusieurs autres occasions, je crois qu’on me permettra d’ajouter encore quelques mots sur cet objet, quoique ce ne soit pas ici le lieu d’en traiter.

Pour présenter la méthode dont il s’agit de la manière la plus simple et en même temps la plus générale qu’il est possible, considérons une ou plusieurs équations différentielles telles que

${\frac {d^{m}x}{dt^{m}}}=\Xi ,\quad {\frac {d^{n}y}{dt^{n}}}=\Psi ,\ldots ,$

dans lesquelles $\Xi ,\Psi ,\ldots$ soient des fonctions données des variables $t,x,y,\ldots$ et des différences ${\frac {dx}{dt}},\ {\frac {dy}{dt}},\ldots ,{\frac {d^{2}x}{dt^{2}}},\ldots$ des ordres inférieurs à ${\frac {d^{m}x}{dt^{m}}},$${\frac {d^{n}y}{dt^{n}}},\ldots .$

Supposons que l’on connaisse une intégrale quelconque de ces équations, laquelle soit représentée par

$\mathrm {V} =0,$

$\mathrm {V}$ étant une fonction des mêmes variables et de leurs différences, et contenant de plus une constante arbitraire $a.$ On sait que, si l’on différentie cette intégrale en y faisant tout varier excepté $a,$ et qu’on y substitue ensuite pour ${\frac {d^{m}x}{dt^{m-1}}},\ {\frac {d^{n}y}{dt^{n-1}}},\ldots$ leurs valeurs $\Xi dt,\Psi dt,\ldots$ tirées des équations différentielles, on doit avoir une équation identique avec l’équation $\mathrm {V} =0,$ ou du moins qui aura lieu en même temps que celle-ci, indépendamment d’aucune relation entre $a,t,x,y,\ldots ,{\frac {dx}{dt}},\ldots \,;$ de sorte que, l’équation $\mathrm {V} =0$ étant posée, la différence $d\mathrm {V}$ deviendra identiquement nulle après les substitutions dont il s’agit.

Soient maintenant proposées les équations

${\frac {d^{m}x}{dt^{m}}}=\Xi +\Xi ',\quad {\frac {d^{n}y}{dt^{n}}}=\Psi +\Psi ',\ldots ,$

$\Xi ',\Psi ',\ldots$ étant de même des fonctions de $t,x,y,\ldots ,{\frac {dx}{dt}},{\frac {dy}{dt}},\ldots ,$ jusqu’à ${\frac {d^{m-1}x}{d^{m-1}t}},{\frac {d^{n-1}y}{d^{n-1}t}},\ldots \,;$ il est clair que la même équation $\mathrm {V} =0$ pourra satisfaire aussi à ces équations, pourvu que la quantité $a,$ étant supposée variable, soit telle que la différentielle $d\mathrm {V}$ devienne nulle en même temps que $\mathrm {V} ,$ après la substitution des valeurs précédentes de ${\frac {d^{m}x}{dt^{m}}},\ {\frac {d^{n}y}{dt^{n}}},\ldots .$ Or nous venons de voir que la partie de $d\mathrm {V} ,$ qui ne contient point la différence $da,$ devient identiquement nulle avec la fonction $\mathrm {V}$ par la substitution des $\Xi ,\Psi ,\ldots$ à la place de ${\frac {d^{m}x}{dt^{m}}},\ {\frac {d^{n}y}{dt^{n}}},\ldots .$ Donc il n’y aura qu’à

rendre nulle la partie restante de $d\mathrm {V} ,$ c’est-à-dire les termes provenant de la différence de $a,$ et de la substitution des quantités $\Xi ',\Psi ',\ldots$ au lieu de ${\frac {d^{m}x}{dt^{m}}},\ {\frac {d^{n}y}{dt^{n}}},\ldots ,$ ce qui donnera équation

${\frac {d\mathrm {V} }{da}}da+{\frac {d\mathrm {V} }{d{\cfrac {d^{m-1}x}{dt^{m-1}}}}}\Xi 'dt+{\frac {d\mathrm {V} }{d{\cfrac {d^{n-1}y}{dt^{n-1}}}}}\Psi 'dt+\ldots =0.$

C’est l’équation qui servira à déterminer la valeur convenable de la quantité $a$ devenue variable.

Si l’on avait deux équations intégrales

$\mathrm {V} =0,\quad \mathrm {U} =0,$

des mêmes équations différentielles

${\frac {d^{m}x}{d^{m}t}}=\Xi ,\quad {\frac {d^{n}y}{d^{n}t}}=\Psi ,\ldots ,$

et que ces intégrales continssent deux constantes arbitraires $a,b,$ on prouverait de la même manière qu’elles pourraient l’être aussi des équations

${\frac {d^{m}x}{d^{m}t}}=\Xi +\Xi ',\quad {\frac {d^{n}y}{d^{n}t}}=\Psi +\Psi ',\ldots ,$

en y supposant $a$ et $b$ variables, et déterminées par ces équations

${\begin{aligned}{\frac {d\mathrm {V} }{da}}da+{\frac {d\mathrm {V} }{db}}db&+{\frac {d\mathrm {V} }{d{\cfrac {d^{m-1}x}{dt^{m-1}}}}}\Xi 'dt+{\frac {d\mathrm {V} }{d{\cfrac {d^{n-1}y}{dt^{n-1}}}}}\Psi 'dt+\ldots =0,\\{\frac {d\mathrm {U} }{da}}da+{\frac {d\mathrm {U} }{db}}db&+{\frac {d\mathrm {U} }{d{\cfrac {d^{m-1}x}{dt^{m-1}}}}}\Xi 'dt\,+{\frac {d\mathrm {U} }{d{\cfrac {d^{n-1}y}{dt^{n-1}}}}}\Psi 'dt+\ldots =0.\end{aligned}}$

Et ainsi de suite si l’ori avait un plus grand nombre d’intégrales.

26. Au reste, quand on connaît une intégrale qui renferme deux constantes arbitraires, on en peut d’abord déduire une seconde par la seule différentiation, en substituant, s’il est nécessaire, à la place des plus hautes différences des variables, leurs valeurs tirées des équations différentielles données. De la même manière une intégrale qui renferme trois constantes arbitraires fournira, par deux différentiations successives, deux autres intégrales, et ainsi de suite. Donc, si les intégrâtes connues renferment autant de constantes arbitraires qu’il y a d’unités dans la somme des exposants des équations différentielles, ce qui est le cas des intégrales finies et complètes, on trouvera par leur moyen autant d’intégrales différentes qu’il y a d’arbitraires ; et l’on aura par la méthode précédente toutes les équations nécessaires pour déterminer les valeurs de ces constantes devenues variables.

De plus on pourra dans ce cas déterminer les variables finies $x,y,\ldots ,$ ainsi que leurs différences ${\frac {dx}{dt}},{\frac {dy}{dt}},\ldots$ jusqu’à ${\frac {d^{m-1}x}{dt^{m-1}}},\ {\frac {d^{n-1}y}{dt^{n-1}}},\ldots ,$ en fonction de $t$ et des constantes arbitraires $a,b,\ldots \,;$ et il est visible que ces fonctions seront les mêmes, soit que ces arbitraires soient variables ou non ; de sorte que, dans les différentiations de $x,y,$ jusqu’à $d^{m-1}x,d^{n-1}y,\ldots ,$ on pourra toujours regarder et traiter les quantités $a,b,c,\ldots$ comme des constantes invariables.

27. La méthode que j’avais donnée dans les Mémoires de 1775^[6] rentre aussi dans celle que je viens d’exposer, et l’on peut généraliser ainsi l’application que j’en avais faite aux équations linéaires. En effet, lorsque $\Xi ,\Psi ,\ldots$ sont des fonctions linéaires des variables $x,y,\ldots$ et de leurs différences, il est facile de prouver que $\mathrm {V,U} ,\ldots$ seront aussi des fonctions linéaires des mêmes variables et des constantes arbitraires $a,b,\ldots .$ Donc ${\frac {d\mathrm {V} }{da}},{\frac {d\mathrm {V} }{db}},{\frac {d\mathrm {V} }{d{\cfrac {d^{m-1}x}{dt}}}},\ldots$ seront des fonctions de $t.$ Par conséquent si $\Xi ',\Psi ',\ldots$ sont données en $t$ seul, on déterminera facilement les valeurs de $a,b,\ldots$ en $t.$

section seconde.

formules générales pour les variations séculaires des éléments
des planètes.

28. Les formules que nous venons de donner, dans la Section précédente, pour représenter les variations des éléments des Planètes, causées par leur action mutuelle, expriment l’effet total de cette action, et pourraient servir à déterminer toutes les inégalités qu’elle doit produire dans leur mouvement. Mais notre objet est simplement de déterminer les variations séculaires des éléments des Planètes, c’est-à-dire celles qui n’ont aucune période fixe, ou du moins qui en ont de très-longues et indépendantes du retour des Planètes aux mêmes points de leurs orbites. Ces variations sont nécessairement renfermées dans les formules trouvées, et, pour les démêler, il n’y aura qu’à développer ces formules et les débarrasser ensuite de tout ce qu’elles peuvent renfermer de périodique. Or la petitesse des excentricités et des inclinaisons des Planètes fait qu’on peut exprimer leurs coordonnées par des séries très-convergentes de sinus et cosinus d’angles proportionnels au temps. Il faudra donc faire ces substitutions à la place de $x,y,z,x',y',z',\ldots ,$ et rejeter ensuite tous les termes qui se trouveront contenir des sinus et des cosinus. Ainsi il faut commencer par chercher les valeurs convenables de $x,y,z$ en $t.$

29. Nous avons déjà remarqué plus haut que, pour faciliter cette recherche, il est à propos d’employer les substitutions

$x=r\cos q,\quad y=r\sin q,$

$r$ étant le rayon de l’orbite projetée sur le plan des $x,y,$ et $q$ l’angle de ce rayon avec l’axe des $x$.

Ces substitutions ont d’ailleurs l’avantage d’être conformes aux usages astronomiques, puisqu’en prenant le plan des $x$ et $y$ pour celui de l’écliptique, et supposant l’axe des $x$ dirigé vers le premier point d’Aries, $r$ sera la distance accourcie de la Planète au Soleil, et $q$ sa longitude héliocentrique.

En mettant ces expressions de $x$ et $y$ dans l’équation

$\mathrm {P} x-\mathrm {Q} y+\mathrm {R} z=0$

du n^o 4, on en tire

$z={\frac {r}{\mathrm {R} }}(\mathrm {Q} \sin q-\mathrm {P} \cos q).$

“De là on aura

$\rho ={\sqrt {x^{2}+y^{2}+z^{2}}}={\frac {r}{\mathrm {R} }}{\sqrt {\mathrm {R} ^{2}+(\mathrm {Q} \sin q-\mathrm {P} \cos q)^{2}}}\,;$

“et l’équation

$g\rho =\mathrm {N} x+\mathrm {M} y+\mathrm {L} z+\Pi ^{2}$

du n^o 7 donnera, en faisant

$\mathrm {B=RN-PL,\quad C=-RM-QL} ,$

comme dans le n^o 8,

$r={\frac {\mathrm {R} \Pi ^{2}}{g{\sqrt {\mathrm {R} ^{2}+(\mathrm {Q} \sin q-\mathrm {P} \cos q)^{2}}}+\mathrm {C} \sin q-\mathrm {B} \cos q}}.$

Et comme les mêmes équations ont lieu aussi en y faisant varier simplement $x,y,z,$ et regardant les autres quantités comme constantes (10), il s’ensuit que, pour avoir les valeurs des différences $dz$ et $dr,$ il suffira de faire varier dans les formules précédentes les quantités $z,r,q,$ en prenant $\mathrm {P,Q,R,\Pi ,B,C}$ pour constantes.

Enfin l’équation

$xdy-ydx=\mathrm {R} dt$

du n^o 3 donnera

$r^{2}dq=\mathrm {R} dt,$

de sorte qu’en substituant pour $r^{2}$ sa valeur en $q$ on aura

${\frac {\mathrm {R} \Pi ^{4}dq}{\left[g{\sqrt {\mathrm {R} ^{2}+(\mathrm {Q} \sin q-\mathrm {P} \cos q)^{2}}}+\mathrm {C} \sin q-\mathrm {B} \cos q\right]^{2}}}=dt.$

Par cette équation oh déterminera donc $q$ en $t\,;$ ensuite on aura, par la substitution de cette valeur de $q,$ celles de $x,y,z$ en $t.$

30. Lorsque l’inclinaison et l’excentricité sont l’une et l’autre fort petites, comme cela a lieu dans notre système planétaire, les quantités $\mathrm {P} \,;$ et $\mathrm {Q}$ sont nécessairement toujours très-petites vis-à-vis de $\mathrm {R} ,$ et les quantités $\mathrm {L,M,N}$ le sont aussi (4, 9), de sorte que $\mathrm {B}$ et $\mathrm {C}$ seront pareillement très-petites par rapport à $\mathrm {R} .$ On pourra donc dans ce cas développer la fraction

${\frac {1}{\left[g{\sqrt {\mathrm {R} ^{2}+(\mathrm {Q} \sin q-\mathrm {P} \cos q)^{2}}}+\mathrm {C} \sin q-\mathrm {B} \cos q\right]^{2}}}$

en une suite fort convergente, laquelle, étant ordonnée relativement aux sinus et cosinus de $q$ et de ses multiples, sera de la forme

$\alpha \left(1+\beta \sin q+\gamma \cos q+\delta \sin 2q+\varepsilon \cos 2q+\ldots \right),$

$\alpha$ étant une quantité finie ; $\beta ,\gamma$ étant des quantités très-petites du premier ordre ; $\delta ,\varepsilon$ étant très-petites du second ordre, et ainsi de suite.

Par ce moyen, l’équation précédente deviendra

$\mathrm {R} \Pi ^{4}\alpha \left(1+\beta \sin q+\gamma \cos q+\delta \sin 2q+\varepsilon \cos 2q+\ldots \right)dq=dt,$

ou bien en divisant par $\mathrm {R} \Pi ^{4}\alpha ,$

$dq+\beta \sin qdq+\gamma \cos qdq+\delta \sin 2qdq+\varepsilon \cos 2qdq+\ldots ={\frac {dt}{\mathrm {R} \Pi ^{4}\alpha }}.$

Le premier membre de cette équation est intégrable exactement lorsque $\beta ,\gamma ,\ldots$ sont des quantités constantes ; mais, lorsque ces quantités sont variables, il faut avoir recourus aux séries ; et l’on trouve ; en employant l’opération connue des intégrations par parties, et regardant $\beta ,\gamma ,\ldots$ comme des fonctions de $q.$

${\begin{aligned}&\int \beta \sin qdq=-\beta \cos q+\int {\frac {d\beta }{dq}}\cos qdq,\\&\int {\frac {d\beta }{dq}}\cos qdq={\frac {d\beta }{dq}}\sin q-\int {\frac {d^{2}\beta }{dq^{2}}}\sin qdq,\\&\int {\frac {d^{2}\beta }{dq^{2}}}\sin qdq=-{\frac {d^{2}\beta }{dq^{2}}}\cos q+\int {\frac {d^{3}\beta }{dq^{3}}}\cos qdq,\end{aligned}}$

et ainsi de suite,

On aura de même

${\begin{aligned}&\int \gamma \cos qdq=\gamma \sin q-\int {\frac {d\gamma }{dq}}\sin qdq,\\&\int {\frac {d\gamma }{dq}}\sin qdq=-{\frac {d\gamma }{dq}}\cos q+\int {\frac {d^{2}\gamma }{dq^{2}}}\cos qdq,\\&\ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots .\end{aligned}}$

Donc, si l’on fait ces substitutions successives, et qu’on suppose, pour abréger,

${\begin{aligned}(\beta )=&\beta -{\frac {d\gamma }{dq}}-{\frac {d^{2}\beta }{dq^{2}}}+{\frac {d^{3}\gamma }{dq^{3}}}+{\frac {d^{4}\beta }{dq^{4}}}-\ldots ,\\(\gamma )=&\gamma +{\frac {d\beta }{dq}}-{\frac {d^{2}\gamma }{dq^{2}}}-{\frac {d^{3}\beta }{dq^{3}}}+{\frac {d^{4}\gamma }{dq^{4}}}+\ldots ,\end{aligned}}$

on aura

$\int (\beta \sin qdq+\gamma \cos qdq)=-(\beta )\cos q+(\gamma )\sin q.$

Et, supposant pareillement.

${\begin{aligned}(\delta )=&\delta -{\frac {1}{2}}{\frac {d\varepsilon }{dq}}-{\frac {1}{4}}{\frac {d^{2}\delta }{dq^{2}}}+{\frac {1}{8}}{\frac {d^{3}\varepsilon }{dq^{3}}}+{\frac {1}{16}}{\frac {d^{4}\delta }{dq^{4}}}-\ldots ,\\(\varepsilon )=&\varepsilon +{\frac {1}{2}}{\frac {d\delta }{dq}}-{\frac {1}{4}}{\frac {d^{2}\varepsilon }{dq^{2}}}-{\frac {1}{8}}{\frac {d^{3}\delta }{dq^{3}}}+{\frac {1}{16}}{\frac {d^{4}\varepsilon }{dq^{4}}}+\ldots ,\end{aligned}}$

on trouvera

$\int (\delta \sin 2qdq+\varepsilon \cos 2qdq)=-{\frac {1}{2}}(\delta )\cos q+{\frac {1}{2}}(\varepsilon )\sin q,$

et ainsi de suite.

À l’égard du second membre de l’équation, il est évidemment intégrable en y regardant $\mathrm {R} ,\Pi ,\alpha$ comme des fonctions de $t.$

Soit, pour plus de simplicité,

$dp={\frac {dt}{\mathrm {R} \Pi ^{4}\alpha }},$

et l’intégrale de l’équation en question sera

$q-(\beta )\cos q+(\gamma )\sin q-{\frac {1}{2}}(\delta )\cos 2q+{\frac {1}{2}}(\varepsilon )\sin 2q+\ldots =p,$

de laquelle, puisque $(\beta ),(\gamma )$ sont supposées très-petites du premier ordre, $(\delta ),(\varepsilon )$ très-petites du second ordre, et ainsi de suite, il est facile de tirer la valeur de $q$ en $p,$ exprimée par une suite fort convergente.

31. En général, si l’on a l’équation

$q+f(q)=p,$

$f(q)$ dénotant une fonction quelconque de $q,$ on aura, par le Théorème que j’ai donné ailleurs,

$q=p-f(p)+{\frac {1}{2}}{\frac {d\left[f(p)\right]^{2}}{dp}}-{\frac {1}{2.3}}{\frac {d^{2}\left[f(p)\right]^{3}}{dp^{2}}}+\ldots ,$

et même, en dénotant par $\varphi$ une autre fonction quelconque et faisant $\varphi '(p)={\frac {d\varphi (p)}{dp}},$

$\varphi (q)=\varphi (p)-f(p)\varphi '(p)+{\frac {1}{2}}{\frac {d\left[\left[f(p)\right]^{2}\varphi '(p)\right]}{dp}}-\ldots .$

Ainsi, dans notre cas, il n’y aura qu’à faire

$f(q)=-(\beta )\cos q+(\gamma )\sin q-\ldots ,$

et par conséquent

$f(p)=-(\beta )\cos p+(\gamma )\sin p-{\frac {1}{2}}(\delta )\cos 2p+{\frac {1}{2}}(\varepsilon )\sin 2p-\ldots ,$

et exécuter ensuite relativement à $p$ les différentiations indiquées.

On trouvera de cette manière en supposant

$q=p-(\mathrm {B} )\cos p-(\mathrm {C} )\sin p+(\mathrm {D} )\cos 2p-(\mathrm {E} )\sin 2p+\ldots ,$

en supposant

${\begin{aligned}(\mathrm {B} )=&(\beta )-{\frac {1}{4}}(\beta )(\varepsilon )+{\frac {1}{4}}(\gamma )(\delta )+\ldots ,\\(\mathrm {C} )=&(\gamma )+{\frac {1}{4}}(\beta )(\delta )+{\frac {1}{4}}(\gamma )(\varepsilon )+\ldots ,\\(\mathrm {D} )=&{\frac {1}{2}}(\delta )-(\beta )(\gamma )+\ldots ,\\(\mathrm {E} )=&{\frac {1}{2}}(\varepsilon )+{\frac {1}{2}}(\beta )^{2}-{\frac {1}{2}}(\gamma )^{2}+\ldots ,\\\ldots \ldots &\ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots .\end{aligned}}$

32. Dans les orbites non troublées la quantité $p$ est proportionnelle au temps $t,$ parce que les quantités $\mathrm {R} ,\Pi ,\alpha$ y sont constantes, en sorte que

$p={\frac {t}{\mathrm {R} \Pi ^{4}\alpha }}.$

Ainsi $p$ est alors la valeur moyenne de $q\,;$ et puisque $q$ est la longitude vraie de la Planète, $p$ en sera la longitude moyenne.

Il n’en est pas de même pour les orbites troublées, ou les quantités $\mathrm {R} ,\Pi ,\alpha$ H, a sont variables ; cependant on peut toujours, par analogie, y regarder la quantité $p$ comme la longitude moyenne ; mais alors le mouvement moyen ne sera plus uniforme, et la vitesse de ce mouvement sejrouvera exprimée par la quantité variable ${\frac {1}{\mathrm {R} \Pi ^{4}\alpha }}.$

Si cette quantité ne contenait que des termes proportionnels aux sinus ou cosinus de $t$ et de ses multiples, il est clair que les variations de $p$ qui en proviendraient ne seraient que périodiques ; elles rentreraient par conséquent dans les inégalités périodiques du mouvement des Planètes, inégalités dont nous faisons abstraction dans ces Recherches. Mais si la quantité ${\frac {1}{\mathrm {R} \Pi ^{4}\alpha }}$ renferme des termes qui croissent en même temps que t, ou qui aient une période très-longue, ces termes donneront des variations séculaires dans le mouvement moyen, et la détermination de ces variations est un des points les plus importants de la Théorie que nous traitons. Il est donc nécessaire de déterminer rigoureusement la loi de la variation de la quantité dont il s’agit, et pour cela il faut connaître la valeur de la quantité qui représente le terme tout constant de la fraction

${\frac {1}{\left[g{\sqrt {\mathrm {R} ^{2}+(\mathrm {Q} \sin q-\mathrm {P} \cos q)^{2}}}+\mathrm {C} \sin q-\mathrm {B} \cos q\right]^{2}}},$

développée suivant les sinus et cosinus des multiples de $q\,;$ c’est de quoi nous allons nous occuper.

33. Commençons par faire disparaître le radical du dénominateur, en multipliant le haut et le bas de la fraction par la quantité

$\left[g{\sqrt {\mathrm {R} ^{2}+(\mathrm {Q} \sin q-\mathrm {P} \cos q)^{2}}}-\mathrm {C} \sin q+\mathrm {B} \cos q\right]^{2},$

on aura cette transformée

${\frac {\left[g{\sqrt {\mathrm {R} ^{2}+(\mathrm {Q} \sin q-\mathrm {P} \cos q)^{2}}}-\mathrm {C} \sin q+\mathrm {B} \cos q\right]^{2}}{\left[g^{2}\mathrm {R} ^{2}+g^{2}(\mathrm {Q} \sin q-\mathrm {P} \cos q)^{2}-(\mathrm {C} \sin q-\mathrm {B} \cos q)^{2}\right]^{2}}}\,;$

laquelle se réduit à cette forme

${\frac {a+b\cos 2q-c\sin 2q+2g(\mathrm {B} \cos q-\mathrm {C} \sin q)\mathrm {V} }{(h+m\cos 2q-n\sin 2q)^{2}}},$

en faisant, pour abréger,

${\begin{aligned}a=&g^{2}\mathrm {\left(R^{2}+{\frac {Q^{2}+P^{2}}{2}}\right)+{\frac {C^{2}+B^{2}}{2}}} ,\\b=&g^{2}\mathrm {{\frac {P^{2}-Q^{2}}{2}}+{\frac {B^{2}-C^{2}}{2}}} ,\\c=&g^{2}\mathrm {PQ+BC} ,\\h=&g^{2}\mathrm {\left(R^{2}+{\frac {Q^{2}+P^{2}}{2}}\right)-{\frac {C^{2}+B^{2}}{2}}} ,\\m=&g^{2}\mathrm {{\frac {P^{2}-Q^{2}}{2}}-{\frac {B^{2}-C^{2}}{2}}} ,\\n=&g^{2}\mathrm {PQ-BC} ,\\\mathrm {V} =&{\sqrt {\mathrm {R} ^{2}+(\mathrm {Q} \sin q-\mathrm {P} \cos q)^{2}}}.\end{aligned}}$

À considérer cette formule, il est facile de voir que la partie qui a pour numérateur $a+b\cos 2q-c\sin 2q$ ne donnera par le développementque des termes proportionnels à des sinus ou cosinus de multiples pairs de $q,$ et que l’autre partie dont le numérateur est $2g(\mathrm {B} \cos q-\mathrm {C} \sin q)\mathrm {V}$ donnera seulement des termes proportionnels aux sinus et cosinus des multiples impairs de $q.$ De sorte qu’on aura (30)

${\begin{aligned}{\frac {a+b\cos 2q-c\sin 2q}{(h+m\cos 2q-n\sin 2q)^{2}}}=&\alpha \left(1+\delta \sin 2q+\varepsilon \cos 2q+\ldots \right),\\{\frac {2g(\mathrm {B} \cos q-\mathrm {C} \sin q)\mathrm {V} }{(h+m\cos 2q-n\sin 2q)^{2}}}=&\alpha \left(\beta \sin q+\gamma \cos q+\zeta \sin 3q+\eta \cos 3q+\ldots \right).\end{aligned}}$

Ainsi la question se réduit à trouver le terme tout constant $\alpha$ de la

fraction rationnelle

${\frac {a+b\cos 2q-c\sin 2q}{(h+m\cos 2q-n\sin 2q)^{2}}},$

développée suivant les sinus et cosinus des multiples de $2q.$

Or, si au lieu de cette fraction on considère celle-ci plus simple

${\frac {a+b\cos 2q-c\sin 2q}{h+m\cos 2q-n\sin 2q}},$

et qu’on la développe en une série de la forme

$\mathrm {A} +\mathrm {D} \sin 2q+\mathrm {E} \cos 2q+\mathrm {F} \sin 4q+\mathrm {G} \cos 4q+\ldots ,$

il est clair qu’en faisant varier de part et d’autre la quantité $h,$ et divisant par $dh,$ on aura

$-{\frac {a+b\cos 2q-c\sin 2q}{(h+m\cos 2q-n\sin 2q)^{2}}}={\frac {d\mathrm {A} }{dh}}+{\frac {d\mathrm {D} }{dh}}\sin 2q+{\frac {d\mathrm {E} }{dh}}\cos 2q+{\frac {d\mathrm {F} }{dh}}\sin 4q+\ldots .$

De sorte qu’on aura par la comparaison des termes

$\alpha =-{\frac {d\mathrm {A} }{dh}},\quad \alpha \delta =-{\frac {d\mathrm {D} }{dh}},\quad \alpha \varepsilon =-{\frac {d\mathrm {E} }{dh}},\ldots .$

Il ne s’agit donc que de développer la dernière fraction ; c’est ce qu’on peut faire par différentes méthodes ; mais aucune ne me paraît plus simple que celle que je vais exposer, et qui peut d’ailleurs être utile aussi dans d’autres occasions.

34. Je fais pour plus de simplicité $2q=u,$ et substituant dans la fraction proposée, à la place de sinu et cosu, leurs valeurs en exponentielles imaginaires, je la réduis à cette forme

${\frac {2a+\left(b+c{\sqrt {-1}}\right)e^{u{\sqrt {-1}}}+\left(b-c{\sqrt {-1}}\right)e^{-u{\sqrt {-1}}}}{2h+\left(m+n{\sqrt {-1}}\right)e^{u{\sqrt {-1}}}+\left(m-n{\sqrt {-1}}\right)e^{-u{\sqrt {-1}}}}}.$

Cette fraction peut se partager en ces deux-ci

${\frac {\lambda +\left(\mu +\nu {\sqrt {-1}}\right)e^{u{\sqrt {-1}}}}{\varpi +\left(\rho +\sigma {\sqrt {-1}}\right)e^{u{\sqrt {-1}}}}}+{\frac {\lambda +\left(\mu -\nu {\sqrt {-1}}\right)e^{-u{\sqrt {-1}}}}{\varpi +\left(\rho -\sigma {\sqrt {-1}}\right)e^{-u{\sqrt {-1}}}}}\,;$

car en multipliant en croix et comparant les termes, on aura ces six équations

${\mathrm{\varpi <!--\; \varpi \; -->}}^2+{\mathrm{\rho <!--\; \rho \; -->}}^2+$