Trajectoires numériques en force centrale

La stratégie. Une jolie courbe ne prouve rien : il faut un critère de justesse. La physique en fournit deux, gratuits : la constante des aires et l'énergie mécanique, que la vraie trajectoire conserve. On calcule d'abord ce qu'on attend, puis on mesure ce que le schéma conserve.

1) Prévisions. En $r_0$ : $v_c=\sqrt{\mathcal{G}M/r_0}=7,55\mathrm{km}{\mathrm{s}}^{-1}$ et $v_{\mathrm{lib}}=\sqrt{2}{v}_c=10,68\mathrm{km}{\mathrm{s}}^{-1}$. On a $v_c<v_0<v_{\mathrm{lib}}$ : l'énergie mécanique est négative, l'état est lié, et la trajectoire est une ellipse dont $r_0$ est le périgée (la vitesse y est orthoradiale et supérieure à la vitesse circulaire, donc le satellite s'éloigne ensuite). Par unité de masse : $$\frac{{\mathcal{E}}_{\mathrm{m}}}{m}=\frac{v_0^2}{2}-\frac{\mathcal{G}M}{r_0}=-1,65\times 10^7\mathrm{J}\mathrm{k}{\mathrm{g}}^{-1},a=-\frac{\mathcal{G}M}{2{\mathcal{E}}_{\mathrm{m}}/m}=1,209\times 10^7\mathrm{m}.$$ La somme des distances au périgée et à l'apogée vaut $2a$ (on le démontrera à l'exercice 5 par les constantes du mouvement), d'où $r_{\max }=2a-r_0=1,72\times 10^7\mathrm{m}$, et la troisième loi de Kepler, admise pour l'ellipse, donne $T=2\pi \sqrt{a^3/(\mathcal{G}M)}=3,67\mathrm{h}$.

2) Ce que conserve chaque schéma. Pour le schéma semi-implicite : $${\overrightarrow{r}}_{n+1}\wedge {\overrightarrow{v}}_{n+1}=({\overrightarrow{r}}_n+{\overrightarrow{v}}_{n+1}\Delta t)\wedge {\overrightarrow{v}}_{n+1}={\overrightarrow{r}}_n\wedge {\overrightarrow{v}}_{n+1},$$ car ${\overrightarrow{v}}_{n+1}\wedge {\overrightarrow{v}}_{n+1}=\overrightarrow{0}$. Puis $${\overrightarrow{r}}_n\wedge {\overrightarrow{v}}_{n+1}={\overrightarrow{r}}_n\wedge {\overrightarrow{v}}_n+\Delta t{\overrightarrow{r}}_n\wedge {\overrightarrow{a}}_n={\overrightarrow{r}}_n\wedge {\overrightarrow{v}}_n,$$ car ${\overrightarrow{a}}_n$ est colinéaire à ${\overrightarrow{r}}_n$ — c'est exactement le caractère central de la force. Le schéma conserve donc la constante des aires à la précision de la machine près, quel que soit $\Delta t$.

Pour le schéma explicite : $${\overrightarrow{r}}_{n+1}\wedge {\overrightarrow{v}}_{n+1}=({\overrightarrow{r}}_n+{\overrightarrow{v}}_n\Delta t)\wedge ({\overrightarrow{v}}_n+{\overrightarrow{a}}_n\Delta t)={\overrightarrow{r}}_n\wedge {\overrightarrow{v}}_n+\Delta t^2{\overrightarrow{v}}_n\wedge {\overrightarrow{a}}_n,$$ et le terme $\Delta t^2{\overrightarrow{v}}_n\wedge {\overrightarrow{a}}_n$ n'est pas nul en général : $C$ dérive, d'autant moins que le pas est petit.

3) Le script et son résultat.

import numpy as np

GM, r0, v0 = 3.99e14, 7.0e6, 9.0e3   # SI ; v0 orthoradiale en r0

def orbite(dt, duree, schema, p=2.0, v=v0):
    x, y, vx, vy = r0, 0.0, 0.0, v
    S = [(x, y, vx, vy)]
    for _ in range(int(duree / dt)):
        k = GM * r0**(p - 2) / np.hypot(x, y)**(p + 1)
        ax, ay = -k * x, -k * y
        if schema == "explicite":     # tout avec l'état ancien
            x, y = x + vx*dt, y + vy*dt
            vx, vy = vx + ax*dt, vy + ay*dt
        else:                         # semi-implicite (le cours)
            vx, vy = vx + ax*dt, vy + ay*dt
            x, y = x + vx*dt, y + vy*dt
        S.append((x, y, vx, vy))
    return np.array(S).T

# énergie par unité de masse, constante des aires
Em = lambda x, y, vx, vy: (vx**2 + vy**2)/2 - GM/np.hypot(x, y)
C = lambda x, y, vx, vy: x*vy - y*vx

a = -GM / (2*Em(r0, 0, 0, v0)); T = 2*np.pi*np.sqrt(a**3/GM)
for schema in ("explicite", "semi-implicite"):
    for dt in (10.0, 1.0):
        s = orbite(dt, T, schema)
        dE = (Em(*s[:, -1]) - Em(*s[:, 0])) / abs(Em(*s[:, 0]))
        dC = (C(*s[:, -1]) - C(*s[:, 0])) / C(*s[:, 0])
        print(f"{schema:14s} dt = {dt:2.0f} s : "
              f"dEm/|Em| = {dE:+.1e} ; dC/C = {dC:+.0e}")

for v in (10.5e3, 10.9e3):            # de part et d'autre de v_lib
    x, y, vx, vy = orbite(10.0, 3e5, "semi", v=v)
    E0 = Em(x[0], y[0], vx[0], vy[0])
    print(f"v0 = {v:.0f} m/s : Em = {E0:+.2e} J/kg ; "
          f"r max sur 3,5 jours = {np.hypot(x, y).max():.2e} m")

x, y, vx, vy = orbite(1.0, 3*T, "semi", p=2.2)   # force en r^-2,2
r, th = np.hypot(x, y), np.unwrap(np.arctan2(y, x))
peri = [k for k in range(100, len(r) - 1) if r[k-1] > r[k] <= r[k+1]]
print("p = 2,2 : angle du premier retour au périgée =",
      np.round(np.degrees(th[peri[0]]), 1), "degrés")

Résultat de l'exécution :

explicite      dt = 10 s : dEm/|Em| = +6.6e-02 ; dC/C = +3e-02
explicite      dt =  1 s : dEm/|Em| = +1.0e-02 ; dC/C = +4e-03
semi-implicite dt = 10 s : dEm/|Em| = -5.2e-05 ; dC/C = -1e-15
semi-implicite dt =  1 s : dEm/|Em| = -5.0e-07 ; dC/C = -1e-14
v0 = 10500 m/s : Em = -1.88e+06 J/kg ; r max sur 3,5 jours = 2.06e+08 m
v0 = 10900 m/s : Em = +2.40e+06 J/kg ; r max sur 3,5 jours = 8.29e+08 m
p = 2,2 : angle du premier retour au périgée = 406.8 degrés

Commentaire. Le schéma explicite fabrique de l'énergie : en une seule période, ${\mathcal{E}}_{\mathrm{m}}$ augmente de $6,6\%$ avec un pas de $10\mathrm{s}$, et encore de $1,0\%$ avec un pas de $1\mathrm{s}$. L'orbite calculée s'élargit à chaque tour : c'est une spirale, pas une ellipse. Le schéma semi-implicite conserve $C$ à $10^{-14}$ près — l'arrondi de la machine, comme l'annonçait la question 2) — et son erreur d'énergie ($5\times 10^{-5}$ puis $5\times 10^{-7}$) est divisée par cent quand le pas est divisé par dix et ne s'accumule pas d'un tour à l'autre, pour le même coût de calcul.

4) Le seuil. À $10,5\mathrm{km}{\mathrm{s}}^{-1}$, juste sous $v_{\mathrm{lib}}=10,68\mathrm{km}{\mathrm{s}}^{-1}$, l'énergie est négative et le satellite plafonne à $r_{\max }=2,06\times 10^8\mathrm{m}$, exactement $2a-r_0$ : ellipse très allongée, de période environ quatre jours. À $10,9\mathrm{km}{\mathrm{s}}^{-1}$, l'énergie est positive et la distance croît sans cesse (au-delà de l'orbite de la Lune après trois jours et demi) : état de diffusion. La transition a lieu à ${\mathcal{E}}_{\mathrm{m}}=0$, c'est-à-dire à $v_0=v_{\mathrm{lib}}$, comme l'annonce la lecture de l'énergie potentielle effective.

5) Une loi de force non newtonienne. Le premier retour au point le plus proche se fait après $406,8^{\circ }$ au lieu de $360^{\circ }$ : le périgée avance de $46,8^{\circ }$ par oscillation radiale, et l'orbite est une rosace qui ne se referme pas. Avec la loi en $1/r^2$, l'avance est nulle à la précision du calcul.

Prévision par linéarisation. Pour une force $-k/r^p$, l'énergie potentielle effective par unité de masse est $\frac{C^2}{2r^2}+\frac{k}{(p-1)r^{p-1}}$. Sa dérivée seconde au minimum $r_c$ (où $C^2/r_c^3=k/r_c^p$) vaut $\frac{3C^2}{r_c^4}-\frac{pk}{r_c^{p+1}}=(3-p)\frac{C^2}{r_c^4}=(3-p){\dot{\theta }}_c^2$. Donc ${\omega }_r=\sqrt{3-p}{\dot{\theta }}_c$, et l'angle entre deux péricentres vaut $\frac{360^{\circ }}{\sqrt{3-p}}$. Pour $p=2$, on retrouve $360^{\circ }$ : orbites fermées. Pour $p=2,2$ : $360^{\circ }/\sqrt{0,8}=402,5^{\circ }$, soit une avance de $42,5^{\circ }$. La simulation donne un peu plus ($46,8^{\circ }$) parce que l'orbite étudiée est assez excentrique ; une orbite presque circulaire lancée à $1,02v_c$ donne $43,3^{\circ }$, conforme à la prévision.

Ce que l'exercice installe. Trois réflexes de la capacité numérique. Calculer avant de simuler, pour savoir ce qu'on doit voir. Tester un schéma sur les constantes du mouvement : c'est le seul contrôle qui ne demande pas la solution exacte. Et préférer un schéma qui respecte la structure du problème — ici la centralité de la force, qui garantit la conservation exacte de $C$. Le résultat de la question 5) est un fait profond : parmi les lois de force en puissance de $r$, seules la loi newtonienne et la force de rappel élastique donnent des orbites liées toujours fermées (théorème de Bertrand) ; l'avance mesurée du périhélie de Mercure a été l'un des premiers indices que la gravitation n'est pas exactement newtonienne.