Résolution d'équations différentielles normales avec Python ~ Gravitation universelle

introduction

Utilisez Python pour résoudre le mouvement gravitationnel universel.

Equation de mouvement (3 corps)

Les équations de mouvement des trois objets qui se déplacent sous l'attraction universelle sont présentées ci-dessous. Les objets se distinguent par les indices s (soleil), e (terre) et l (luna).

Masse

m_s \\
m_e \\
m_l \\

Vecteur de position

\vec{r_s} = (r_{sx}, r_{sy}, r_{sz}) \\
\vec{r_e} = (r_{ex}, r_{ey}, r_{ez}) \\
\vec{r_l} = (r_{lx}, r_{ly}, r_{lz}) \\

Vecteur de position relative

\vec{r_{se}} = \vec{r_e} - \vec{r_s} = - \vec{r_{es}} \\
\vec{r_{sl}} = \vec{r_l} - \vec{r_s} = - \vec{r_{ls}} \\
\vec{r_{el}} = \vec{r_l} - \vec{r_e} = - \vec{r_{le}} \\
r_{se} = | \vec{r_{se}} | = r_{es} \\
r_{sl} = | \vec{r_{sl}} | = r_{ls} \\
r_{el} = | \vec{r_{el}} | = r_{le} \\

Vecteur de vitesse

\vec{v_s} = (v_{sx}, v_{sy}, v_{sz}) = \dot{\vec{r_s}} \\
\vec{v_e} = (v_{ex}, v_{ey}, v_{ez}) = \dot{\vec{r_e}} \\
\vec{v_l} = (v_{lx}, v_{ly}, v_{lz}) = \dot{\vec{r_l}} \\

Équation du mouvement

\vec{F_s} = (F_{sx}, F_{sy}, F_{sz}) = m_s \ddot{\vec{r_s}} = m_s \dot{\vec{v_s}} \\
\vec{F_e} = (F_{ex}, F_{ey}, F_{ez}) = m_e \ddot{\vec{r_e}} = m_e \dot{\vec{v_e}} \\
\vec{F_l} = (F_{lx}, F_{ly}, F_{lz}) = m_l \ddot{\vec{r_l}} = m_l \dot{\vec{v_l}} \\

Gravitation universelle

\vec{F_s} = \vec{F_{se}} + \vec{F_{sl}} \\
\vec{F_e} = \vec{F_{es}} + \vec{F_{el}} \\
\vec{F_l} = \vec{F_{ls}} + \vec{F_{le}} \\
\vec{F_{se}} = - \vec{F_{es}} = G * m_s * m_e \frac{\vec{r_{se}}}{r_{se}^3} \\
\vec{F_{sl}} = - \vec{F_{ls}} = G * m_s * m_l \frac{\vec{r_{sl}}}{r_{sl}^3} \\
\vec{F_{el}} = - \vec{F_{le}} = G * m_e * m_l \frac{\vec{r_{el}}}{r_{el}^3} \\

Equation de mouvement de 3 objets sous attraction universelle

m_s * \ddot{\vec{r_s}} = G * m_s * m_e \frac{\vec{r_{se}}}{r_{se}^3} + G * m_s * m_l \frac{\vec{r_{sl}}}{r_{sl}^3} \\
\dot{\vec{v_s}} = \ddot{\vec{r_s}} = G * m_e \frac{\vec{r_{se}}}{r_{se}^3} + G * m_l \frac{\vec{r_{sl}}}{r_{sl}^3} \\

m_e * \ddot{\vec{r_e}} = G * m_e * m_s \frac{\vec{r_{es}}}{r_{es}^3} + G * m_e * m_l \frac{\vec{r_{el}}}{r_{el}^3} \\
\dot{\vec{v_e}} = \ddot{\vec{r_e}} = G * m_s \frac{\vec{r_{es}}}{r_{es}^3} + G * m_l \frac{\vec{r_{el}}}{r_{el}^3} \\
\dot{\vec{v_e}} = \ddot{\vec{r_e}} = - G * m_s \frac{\vec{r_{se}}}{r_{se}^3} + G * m_l \frac{\vec{r_{el}}}{r_{el}^3} \\

m_l * \ddot{\vec{r_l}} = G * m_l * m_s \frac{\vec{r_{ls}}}{r_{ls}^3} + G * m_l * m_e \frac{\vec{r_{le}}}{r_{le}^3} \\
\dot{\vec{v_l}} = \ddot{\vec{r_l}} = G * m_s \frac{\vec{r_{ls}}}{r_{ls}^3} + G * m_e \frac{\vec{r_{le}}}{r_{le}^3} \\
\dot{\vec{v_l}} = \ddot{\vec{r_l}} = - G * m_s \frac{\vec{r_{sl}}}{r_{sl}^3} - G * m_e \frac{\vec{r_{el}}}{r_{el}^3} \\

Format de scipy.integrate.ode

L'équation cinétique de la gravitation universelle est exprimée dans un format géré par scipy.integrate.ode.

\vec{y} = (\vec{r_{s}}, \vec{r_{e}}, \vec{r_{l}}, \vec{v_{s}}, \vec{v_{e}}, \vec{v_{l}}) \\
\dot{\vec{y}} = (\dot{\vec{r_{s}}}, \dot{\vec{r_{e}}}, \dot{\vec{r_{l}}}, \dot{\vec{v_{s}}}, \dot{\vec{v_{e}}}, \dot{\vec{v_{l}}}) \\
\dot{\vec{y}} = (\vec{v_{s}}, \vec{v_{e}}, \vec{v_{l}}, \dot{\vec{v_{s}}}, \dot{\vec{v_{e}}}, \dot{\vec{v_{l}}}) \\

\dot{\vec{v_{s}}} = G * m_e \frac{\vec{r_{se}}}{r_{se}^3} + G * m_l \frac{\vec{r_{sl}}}{r_{sl}^3} \\
\dot{\vec{v_{e}}} = - G * m_s \frac{\vec{r_{se}}}{r_{se}^3} + G * m_l \frac{\vec{r_{el}}}{r_{el}^3} \\
\dot{\vec{v_{l}}} = - G * m_s \frac{\vec{r_{sl}}}{r_{sl}^3} - G * m_e \frac{\vec{r_{el}}}{r_{el}^3} \\

Valeur spécifique

Les masses, distances, vitesses et constantes physiques du soleil, de la terre et de la lune sont:

Python Program

test.py

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from scipy.integrate import ode

G = 6.67E-11
ms = 1.99E+30
me = 5.97E+24
ml = 7.35E+22
Gms = G*ms
Gme = G*me
Gml = G*ml
re0 = 1.50E+11
rl0 = 3.84E+08
ve0 = 2.98E+04
vl0 = 1.02E+03

# initial condition
# [0]:rsx,  [1]:rsy,  [2]:rsz,  [3]:rex,  [4]:rey,  [5]:rez, [6]:rlx,  [7]:rly,  [8]:rlz,  
# [9]:vsx, [10]:vsy, [11]:vsz, [12]:vex, [13]:vey, [14]:vez,[15]:vlx, [16]:vly, [17]:vlz,  
y0 = np.array([
  0,    # rsx
  0,    # rsy
  0,    # rsz
  re0,  # rex
  0,    # rey
  0,    # rez
  re0,  # rlx
  0,    # rly
  rl0,  # rlz
  0,    # vsx
  0,    # vsy
  0,    # vsz
  0,    # vex
  ve0,  # vey
  0,    # vez
  vl0,  # vlx
  ve0,  # vly
  0,    # vlz
])

# y
# [0]:rsx,  [1]:rsy,  [2]:rsz,  [3]:rex,  [4]:rey,  [5]:rez, [6]:rlx,  [7]:rly,  [8]:rlz,  
# [9]:vsx, [10]:vsy, [11]:vsz, [12]:vex, [13]:vey, [14]:vez,[15]:vlx, [16]:vly, [17]:vlz,  
# return
# [0]:rsx',  [1]:rsy',  [2]:rsz',  [3]:rex',  [4]:rey',  [5]:rez', [6]:rlx',  [7]:rly',  [8]:rlz',  
# [9]:vsx', [10]:vsy', [11]:vsz', [12]:vex', [13]:vey', [14]:vez',[15]:vlx', [16]:vly', [17]:vlz',  
def f(t, y):
    rse = np.array([y[3]-y[0], y[4]-y[1], y[5]-y[2]]) # sun to earth
    nse = np.linalg.norm(rse)
    rsl = np.array([y[6]-y[0], y[7]-y[1], y[8]-y[2]]) # sun to luna
    nsl = np.linalg.norm(rsl)
    rel = np.array([y[6]-y[3], y[7]-y[4], y[8]-y[5]]) # earth to luna
    nel = np.linalg.norm(rel)
    vsdot =   Gme * rse / nse**3 + Gml * rsl / nsl**3
    vedot = - Gms * rse / nse**3 + Gml * rel / nel**3
    vldot = - Gms * rsl / nsl**3 - Gme * rel / nel**3
    return [y[9], y[10], y[11], y[12], y[13], y[14], y[15], y[16], y[17],
            vsdot[0], vsdot[1], vsdot[2],
            vedot[0], vedot[1], vedot[2],
            vldot[0], vldot[1], vldot[2]]

solver = ode(f)
solver.set_integrator(name="dop853")
solver.set_initial_value(y0)

dt = 60 * 60 * 24
tw = dt * 365 * 2
ts = []
rsx = []
rsy = []
rsz = []
rex = []
rey = []
rez = []
rlx = []
rly = []
rlz = []
while solver.t < tw:
    solver.integrate(solver.t+dt)
    m = ms + me + ml;
    # center of mass
    com = [1/m * (ms * solver.y[0] + me * solver.y[3] + ml * solver.y[6]),
           1/m * (ms * solver.y[1] + me * solver.y[4] + ml * solver.y[7]),
           1/m * (ms * solver.y[2] + me * solver.y[5] + ml * solver.y[8])]
    ts += [solver.t]
    rsx += [solver.y[0] - com[0]]
    rsy += [solver.y[1] - com[1]]
    rsz += [solver.y[2] - com[2]]
    rex += [solver.y[3] - com[0]]
    rey += [solver.y[4] - com[1]]
    rez += [solver.y[5] - com[2]]
    rlx += [solver.y[6] - com[0]]
    rly += [solver.y[7] - com[1]]
    rlz += [solver.y[8] - com[2]]

plt.figure(0)
plt.axes(projection='3d')
plt.plot(rsx, rsy, rsz, "r-+")
plt.plot(rex, rey, rez, "b-+")
plt.plot(rlx, rly, rlz, "g-+")
plt.show()

Résultat d'exécution