Lösen normaler Differentialgleichungen mit Python ~ Universal Gravitation

Einführung

Verwenden Sie Python, um die universelle Gravitationsbewegung zu lösen.

Bewegungsgleichung (3 Körper)

Die Bewegungsgleichungen der drei Objekte, die sich unter universeller Anziehungskraft bewegen, sind unten gezeigt. Objekte werden durch die Indizes s (Sonne), e (Erde) und l (Luna) unterschieden.

Masse

m_s \\
m_e \\
m_l \\

Positionsvektor

\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}) \\

Relativer Positionsvektor

\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} \\

Geschwindigkeitsvektor

\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}} \\

Bewegungsgleichung

\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}} \\

Universale Gravitation

\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} \\

Bewegungsgleichung von 3 Objekten, die einer universellen Anziehung unterliegen

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 von scipy.integrate.ode

Die kinetische Gleichung der universellen Gravitation wird in einem Format ausgedrückt, das von scipy.integrate.ode verarbeitet wird.

\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} \\

Spezifischer Wert

Die Massen, Entfernungen, Geschwindigkeiten und physikalischen Konstanten von Sonne, Erde und Mond sind:

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()

Ausführungsergebnis