Programmation scientifique Collection Petit Tech en Python

À propos de la création d'Integrator (Propagator)

Lors de la création d'un intégrateur en Python, la fonction solve_ivp () de scipy est pratique. Dans ce cas, autorisez numpy pour la définition de fonction de ODE. Par exemple, dans le cas du problème de restriction de cercle à trois corps

\begin{align}
    \ddot{x} -2\dot{y} 
        &= x -\dfrac{1-\mu}{r_{1}^3} (x+\mu) + \dfrac{\mu}{r_{2}^3} (1-\mu-x)
    \\
    \ddot{y} +2\dot{x} 
        &= y -\dfrac{1-\mu}{r_{1}^3}y - \dfrac{\mu}{r_{2}^3}y
    \\
    \ddot{z} 
        &= -\dfrac{1-\mu}{r_{1}^3}z - \dfrac{\mu}{r_{2}^3}z
\end{align}

Donc, si vous exprimez cela comme 6 ODE de premier ordre

\begin{align}
    \dfrac{d x}{d t} &= \dot{x}
\\
    \dfrac{d y}{d t} &= \dot{y}
\\
    \dfrac{d z}{d t} &= \dot{z}
\\
    \dfrac{d \dot{x}}{d t} 
        &= \dfrac{d^2 x}{d t^2} = 2\dot{y} +x -\dfrac{1-\mu}{r_{1}^3} (x+\mu) + \dfrac{\mu}{r_{2}^3} (1-\mu-x)
    \\
    \dfrac{d \dot{y}}{d t} 
        &= \dfrac{d^2 y}{d t^2} =-2\dot{x} + y -\dfrac{1-\mu}{r_{1}^3}y - \dfrac{\mu}{r_{2}^3}y
    \\
    \dfrac{d \dot{z}}{d t} 
        &= \dfrac{d^2 z}{d t^2} =-\dfrac{1-\mu}{r_{1}^3}z - \dfrac{\mu}{r_{2}^3}z
\end{align}

Quand c'est une fonction

import numpy as np
from numba import jit

# deifne RHS function in CR3BP
@jit(nopython=True)
def rhs_cr3bp(t,state,mu):
    """Equation of motion in CR3BP, formulated for scipy.integrate.solve=ivp(), compatible with njit
    
    Args:
        t (float): time
        state (np.array): 1D array of Cartesian state, length 6
        mu (float): CR3BP parameter
    Returns:
        (np.array): 1D array of derivative of Cartesian state    
    """

    # unpack positions
    x  = state[0]
    y  = state[1]
    z  = state[2]
    # unpack velocities
    vx = state[3]
    vy = state[4]
    vz = state[5]
    # compute radii to each primary
    r1 = np.sqrt( (x+mu)**2 + y**2 + z**2 )
    r2 = np.sqrt( (x-1+mu)**2 + y**2 + z**2 )
    # setup vector for dX/dt
    deriv = np.zeros((6,))
    # position derivatives
    deriv[0] = vx
    deriv[1] = vy
    deriv[2] = vz
    # velocity derivatives
    deriv[3] = 2*vy + x - ((1-mu)/r1**3)*(mu+x) + (mu/r2**3)*(1-mu-x)
    deriv[4] = -2*vx + y - ((1-mu)/r1**3)*y - (mu/r2**3)*y
    deriv[5] = -((1-mu)/r1**3)*z - (mu/r2**3)*z
    
    return deriv

Sera. Comme dans le terme précédent, dans la plupart des cas, l'expression ODE semble être relativement facile à convertir en numba, c'est pourquoi elle est recommandée. Après cela, vous pouvez appeler la fonction `` résoudre_ivp () ''.

from scipy.integrate import solve_ivp

# integrate until final time tf
tf = 3.05
# initial state
state0 = np.array([ 9.83408400e-01, -9.42453366e-04, 1.27227988e-03, 
                    7.03724138e-01, -1.78296421e+00, 1.13566847e+00])
# cr3bp mass parameter mu
mu = 0.012150585609624
# call solve_ivp
sol = solve_ivp(fun=rhs_cr3bp, t_span=(0,tf), y0=state0, args=(mu,))

Ici, args contient les paramètres requis autres que la matrice d'état ODE (état) et l'heure (t) (dans l'exemple ci-dessus, `` `mu``` est entré. Puisqu'il s'agit d'un tuple, même d'une entrée N'oubliez même pas une virgule (,).