Animation d'une orbite

1. Problème à 2 corps

Nous souhaitons réaliser une animation de la Terre en orbite autour du Soleil. Nous savons depuis Kepler que celle-ci décrit une ellipse dont le Soleil est un foyer et que la Terre accélère quand elle s'approche de son étoile.

La trajectoire peut être décrite en fonction du temps

en coordonnées cartésiennes,
- x(t) = a.cos t
- y(t) = b.sin t

ou en coordonnées polaires.
- r(t) = p / (1 + e . θ(t))
- θ(t) = t

Ces équations sont en accord avec la première loi de Kepler, mais pas avec la seconde. En effet, t augmente de façon linéaire et donc la vitesse de la Terre sur cette orbite est constante.

On pourrait alors remplacer les équations précédentes par celles-ci et tenter de déterminer la fonction f(t).

en cooordonnées cartésiennes,
- x(t) = a.cos f(t)
- y(t) = b.sin f(t)

ou en coordonnées polaires.
- r(t) = p / (1 + e . θ(t))
- θ(t) = f(t)

Pour cela, il faut revenir à la loi de Newton qui dit que la somme des forces est égale à la masse multipliée par l'accélération. L'accélération est la dérivée de la vitesse qui est elle-même dérivée de la trajectoire. On obtient alors :

m . x''(t) = f_x

m . y''(t) = f_y

Où m est la masse de la Terre et (f_x, f_y) les coordonnées du vecteur de la force de gravitation. Cette dernière est égale à une constante multiplié par la masse de la Terre, multipliée par la masse du Soleil, le tout divisé par la distance Terre-Soleil au carré.

Notre équation devient alors :

m . x''(t) = \frac{-G m M x(t)}{(x^2 + y^2)^{2/3}}

m . x''(t) = \frac{-G m M y(t)}{(x^2 + y^2)^{2/3}}

On constate que la masse de la Terre se simplifie. Cela signifie que l'orbite ne dépend que de la masse du Soleil. En effet, si la Terre était coupée en deux, les deux morceaux se suivraient sur la même orbite.

Malheureusement, cela nous mène à une équation différentielle de second ordre qui n'est pas facile à résoudre. On renonce donc à utiliser des équations pour simuler l'orbite de la Terre sur ordinateur. Mais cette petite excursion par la loi de Newton nous donne l'idée d'utiliser les vecteurs pour notre petite animation.

Pour simplifier, on place le Soleil au point (0, 0) et on suppose qu'il est immobile. On va utiliser 6 variables qui vont varier à chaque itération de notre algorithme, qui correspond à une image de l'animation et à une unité de temps t :

x, y : corrdonnées de la Terre à l'instant t.

v_x, v_y : vitesse de la Terre à l'instant t.

a_x, a_y : accéleration de la Terre à l'instant t.

On arrive alors à l'algorithme suivant :

def compute(self, x, y, vx, vy, time, step=1.0):
        t = 0.0
        coords = [(x,y)]
        G = 900.0
        while t <= time:
            r2 = x * x + y * y
            r = math.sqrt(r2)
            r3 = r * r2
            ax = -G * x / r3
            ay = -G * y / r3
            vx += ax * step
            vy += ay * step
            x += vx * step
            y += vy * step
            coords.append((x,y))
            t += step
        return coords

On constate que la seconde loi de Kepler est bien respectée. Mais on voit que la précision n'est pas exceptionnelle. Si vous attendez quelques révolutions, vous verrez que la Terre quitte peu à peu son orbite. En fait, elle reste sur une ellipse, mais celle-ci subit une lente rotation de son grand axe.

Cette imprécision se comprend mieux avec le tracé rouge. On a juste changé la valeur de pas (step) et on on voit que la rotation de l'axe de l'ellipse s'est accéléré. La bonne nouvelle, cependant, c'est que l'ellipse conserve ses caractéristiques (grand et petit diamètres) malgrés l'imprécision.

2. Problème à n corps.

Autant, la résolution mathématique du problème à n corps est ardue, autant il est facile de le simuler sur ordinateur (avec les imprécisions dont on a déjà parlé). En fait, chaque corps subira une accélération due aux n-1 autres. Les corps les plus massifs exerceront une force plus grande.

En voici un exemple.