Calcul de l'orbite d'une planète du système solaire avec Excel

• Données initiales : t = 0, x(0), y(0), v_x(0), v_y(0)
• Choix : le soleil est situé au point (0,0) et ne bouge pas.
• On prend un pas de temps constant : Δt
• calcul 2D
• Algorithme :
  • on calcule l'accélération de la planète au temps t : a = F / m
    • a_x(t) = − G M x(t) / r³(t)
    • a_y(t) = − G M y(t) / r³(t)
  • on calcule la première vitesse au temps Δt/2 :
    • v_x(Δt/2) = v_x(0) + a_x(0) * Δt / 2
    • v_y(Δt/2) = v_y(0) + a_y(0) * Δt / 2
  • on calcule les vitesses suivantes décalées aux temps t+Δt/2 : a = dv/dt ≈ [ v(t+Δt/2) − v(t−Δt/2) ] / Δt
    • v_x(t+Δt/2) = v_x(t−Δt/2) + a_x(t) * Δt
    • v_y(t+Δt/2) = v_y(t−Δt/2) + a_y(t) * Δt
  • on calcule la position suivante au temps t+Δt : v = dr/dt ≈ [ r(t+Δt) − r(t) ] / Δt
    • x(t+Δt) = x(t) + v_x(t+Δt/2) * Δt
    • y(t+Δt) = y(t) + v_y(t+Δt/2) * Δt
• Programmation Excel : fichier Excel complet orbite.xls
  • On réserve quelques couples de cases pour les données : Libellés : A1:A6 et Valeurs : B1:B6
  • ligne 1 : Libellés C1:J1
  • ligne 2 : valeurs de départ : x, y, r, ax, ay au temps t=0, vx(Δt/2), vy(Δt/2)
  • ligne 3 : cases C3:J3 de la ligne 3 à recopier vers le bas

  	A	B	C	D	E	F	G	H	I	J
  1	x(0)	1	t	x(t)	y(t)	r(t)^3	ax(t)	ax(t)	vx(t+dt/2)	vy(t+dt/2)
  2	y(0)	0	0	= $B$1	= $B$2	= racine(d2*d2 + e2*e2)^3	= -$B$6 * d2 / f2	= -$B$6 * e2 / f2	= $B$3 + g2 * $b$5 / 2	= $B$4 + h2 * $b$5 / 2
  3	vx(0)	0	= c2 + $B$5	= d2 + i2 * $B$5	= e2 + j2 * $B$5	= racine(d3*d3 + e3*e3)^3	= -$B$6 * d3 / f3	= -$B$6 * e3 / f3	= i2 + g3 * $b$5	= j2 + h3 * $b$5
  4	vy(0)	1
  5	dt	0.1
  6	GM	1

• affichage de l'orbite :
  • On sélectionne les colonnes (x,y) : d et e
  • On insère une image : on choisit la représentation nuage de points.
  • pour les petites vitesses ( v² / 2 − GM / r < 0 ) ; si la courbe ne représente pas une ellipse qui repasse sur elle-même après un tour, c'est que le pas de temps est trop grand.
• Programmation Python : ( on peut facilement ajouter d'autres planètes et prendre en compte leurs interactions )

  # orbite de la planete
  # donnees
      x0 = 1.
      y0 = 0.
      vx0 = 0.1
      vy0 = 1.3
      dt = 0.1
      GM = 1.
      n = 400
  # tableau des resultats
      x = (n+1) * [0.]
      y = (n+1) * [0.]
  # temps 0 pour x, y, vx, vy
      x[0] = x0
      y[0] = y0
      vx = vx0
      vy = vy0
      # premier demi pas de temps pour passer de v(0) a v(1/2)
      dt_vitesse = dt / 2
      for i in range(n) :
          # r(n) et acceleration(n)
          r2 = x[i] * x[i] + y[i] * y[i]
          r = math.sqrt( r2 )
          ax = - GM * x[i] / (r2 * r)
          ay = - GM * y[i] / (r2 * r)
          # v(n+1/2)
          vx = vx + ax * dt_vitesse
          vy = vy + ay * dt_vitesse
          # r(n+1)
          x[i+1] = x[i] + vx * dt
          y[i+1] = y[i] + vy * dt
          # pour les pas de temps suivants
          dt_vitesse = dt
      # affichage
      orbite0 = pylab.plot(x,y)
      pylab.show()