Chapitre 9 : Résolution d'équations différentielles avec SciPy

Objectifs du chapitre

1. SciPy et le module Integrate

L'outil que nous allons utiliser est odeint (Ordinary Differential Equation Integrator). C'est le standard industriel pour résoudre des systèmes dynamiques.

# Importation standard pour la mécanique

import matplotlib.pyplot as plt
from scipy.integrate import odeint # L'outil magique 

Comment "pense" odeint ?

Pour utiliser odeint, tu dois comprendre que l'ordinateur ne connaît pas la physique. Il ne connaît pas "la gravité" ou "les frottements". Il ne connaît qu'une seule chose : La Pente (la dérivée).

La résolution numérique fonctionne toujours selon le même schéma en 3 étapes :

• L'État actuel (S) : "Je suis ici, à cette vitesse."
• La Dérivée (S′) : "Selon les lois de la physique (Newton), je varie de tant."
• L'Intégration : "Je calcule le point suivant."

La syntaxe

La fonction odeint s'utilise toujours avec 3 arguments principaux.

Le concept du Vecteur d'État (S)

Une équation différentielle mécanique est souvent du 2ème ordre (elle contient une accélération x¨). Or, odeint ne sait résoudre que des équations du 1er ordre (y′=f(y)).

L'astuce consiste à grouper la position et la vitesse dans un seul vecteur appelé Vecteur d'État S.

Pour un mouvement 1D (comme la chute libre), on pose : $$ S = \begin{bmatrix} position \ vitesse \end{bmatrix} $$

On doit alors fournir à Python la dérivée de ce vecteur (S′ ou Sp dans le code) : $$ S' = \frac{dS}{dt} = \begin{bmatrix} \dot{position} \ \dot{vitesse} \end{bmatrix} = \begin{bmatrix} vitesse \ accélération \end{bmatrix} $$

3. Mise en pratique : Résolution d'une équation simple

Prenons l'exemple du boulet de canon ou d'une chute avec frottement.

Supposons l'équation : $\dot{v} = -g - \frac{k}{m}v$.

Voici la structure obligatoire d'un script de résolution :

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# ==========================================
# 1. DÉFINITION DES PARAMÈTRES PHYSIQUES
# ==========================================
g = 9.81  # Accélération de la pesanteur (m/s²)
k = 0.5   # Coefficient de frottement fluide (kg/s)
m = 10.0  # Masse de l'objet (kg)

# ==========================================
# 2. LE MODÈLE PHYSIQUE (FONCTION DÉRIVÉE)
# ==========================================
def modele_chute(vecteur_etat, t):
    """
    Modèle pour un mouvement vertical 1D.
    Entrée : vecteur_etat = [Vitesse actuelle]
    Sortie : [Accélération] (qui est la dérivée de la vitesse)
    """
    
    # A. Lecture de l'état actuel
    v = vecteur_etat[0] # Vitesse instantanée
    
    # B. Application du PFD (Newton)
    # Somme des forces = Poids + Frottement
    # Poids : toujours vers le bas (-g)
    # Frottement : s'oppose à la vitesse (-k/m * v)
    
    # Note : Si v est positif (montée), le terme -k*v est négatif (freinage vers le bas).
    #        Si v est négatif (chute), le terme -k*v devient positif (freinage vers le haut).
    acceleration = -g - (k/m)*v
    
    # C. Retour de la dérivée
    return [acceleration]

# ==========================================
# 3. PRÉPARATION DE LA SIMULATION
# ==========================================

# A. Conditions Initiales
v0 = 50.0  # On lance l'objet vers le haut à 50 m/s
etat_initial = [v0]

# B. Axe temporel
t = np.linspace(0, 10, 100) # Simulation sur 10 secondes

# ==========================================
# 4. RÉSOLUTION NUMÉRIQUE
# ==========================================
solution = odeint(modele_chute, etat_initial, t)

# Extraction des données pour la clarté
# solution[:, 0] contient la première (et unique) variable : la vitesse
vitesse_simulee = solution[:, 0] 

# ==========================================
# 5. VISUALISATION SCIENTIFIQUE
# ==========================================
plt.figure(figsize=(8, 5))

# Tracé de la courbe
plt.plot(t, vitesse_simulee, label="Vitesse de l'objet", color='blue', linewidth=2)

# Ajout d'une ligne de référence à 0 (l'apogée)
plt.axhline(0, color='black', linewidth=1, linestyle='-', label="Apogée (v=0)")

# Habillage "Ingénieur"
plt.title("Chute verticale avec frottements fluides")
plt.xlabel("Temps (s)")
plt.ylabel("Vitesse verticale (m/s)")
plt.grid(True)
plt.legend()

plt.show()

4. Résolution du 2nd Ordre (Exemple du Pendule)

Pour le pendule (et le planeur), on a une position et une vitesse. Le vecteur $S$ a donc deux cases.

L'équation du pendule est : $\ddot{\theta} = -\frac{g}{L}\sin(\theta)$.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# ==========================================
# 1. DÉFINITION DES PARAMÈTRES PHYSIQUES
# ==========================================
g = 9.81  # Accélération pesanteur (m/s²)
L = 1.0   # Longueur du fil (m)

# ==========================================
# 2. LE MODÈLE PHYSIQUE (FONCTION DÉRIVÉE)
# ==========================================
# Cette fonction traduit les équations de Newton pour Python
def modele_pendule(vecteur_etat, t):
    """
    Entrée : 
       vecteur_etat : Liste [Position actuelle, Vitesse actuelle]
       t            : Le temps (non utilisé ici, mais obligatoire pour odeint)
    Sortie :
       [Dérivée de la position, Dérivée de la vitesse]
    """
    
    # A. On déballe le colis (Lecture de l'état actuel)
    theta = vecteur_etat[0]  # Angle actuel (en radians)
    omega = vecteur_etat[1]  # Vitesse angulaire actuelle (en rad/s)
    
    # B. On applique la Physique (Calcul des variations)
    
    # Variation 1 : Comment change la position ?
    # -> La variation de l'angle, c'est la vitesse angulaire.
    d_theta_dt = omega
    
    # Variation 2 : Comment change la vitesse ?
    # -> La variation de la vitesse, c'est l'accélération (Formule du pendule)
    d_omega_dt = -(g / L) * np.sin(theta)
    
    # C. On renvoie les dérivées pour que odeint calcule le pas suivant
    return [d_theta_dt, d_omega_dt]

# ==========================================
# 3. PRÉPARATION DE LA SIMULATION
# ==========================================

# A. Conditions Initiales (Ce qu'on ferait en labo)
angle_depart_degres = 10
vitesse_depart = 0

# Attention : Python calcule TOUJOURS en radians
theta_0 = np.radians(angle_depart_degres) 
omega_0 = vitesse_depart

etat_initial = [theta_0, omega_0]

# B. L'axe temporel (Pendant combien de temps on observe ?)
# 200 photos prises sur 5 secondes
t = np.linspace(0, 5, 200)

# ==========================================
# 4. RÉSOLUTION NUMÉRIQUE (Le Calcul)
# ==========================================
solution = odeint(modele_pendule, etat_initial, t)

# On extrait les résultats pour plus de clarté
positions_rad = solution[:, 0]
vitesses_rad  = solution[:, 1]

# Conversion en degrés pour l'affichage (plus parlant pour l'humain)
positions_deg = np.degrees(positions_rad)

# ==========================================
# 5. VISUALISATION SCIENTIFIQUE (Le Rapport)
# ==========================================

# Création d'une planche avec 2 graphiques (l'un sous l'autre)
# sharex=True permet de zoomer sur les deux en même temps
fig, (ax_pos, ax_vit) = plt.subplots(2, 1, figsize=(8, 8), sharex=True)

# --- Graphique 1 : Position ---
ax_pos.plot(t, positions_deg, color='red', linewidth=2)
ax_pos.set_title("Cinématique du pendule simple (L=1m)")
ax_pos.set_ylabel("Angle (degrés)") # Unité claire
ax_pos.grid(True)

# --- Graphique 2 : Vitesse ---
ax_vit.plot(t, vitesses_rad, color='blue', linestyle='--')
ax_vit.set_ylabel("Vitesse angulaire (rad/s)") # Unité claire
ax_vit.set_xlabel("Temps (s)") # L'abscisse n'est nécessaire qu'en bas
ax_vit.grid(True)

# Ajustement automatique des marges
plt.tight_layout()
plt.show()

5. Systèmes complexes (Exemple : Le Planeur)

Pour un planeur, le mouvement se fait en X et en Y. Tu as donc 4 variables à suivre simultanément : $$S = [x, y, v_x, v_y]$$

Ta fonction Sp devra donc retourner 4 valeurs : $$S' = [v_x, v_y, a_x, a_y]$$

C'est exactement la même logique, mais avec plus de lignes.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# ==========================================
# 1. DÉFINITION DES PARAMÈTRES PHYSIQUES
# ==========================================
g = 9.81       # Gravité (m/s²)
rho = 1.225    # Densité de l'air (kg/m³)
m = 300.0      # Masse du planeur (kg)

# Caractéristiques aérodynamiques
S_aile = 10.0  # Surface des ailes (m²)
Cz = 0.7       # Coefficient de Portance (Lift)
Cx = 0.1       # Coefficient de Traînée (Drag) - incluant la surface frontale

# Pré-calcul des coefficients aérodynamiques constants
k_portance = 0.5 * rho * S_aile * Cz
k_trainee  = 0.5 * rho * S_aile * Cx   # Note : ici on simplifie en utilisant S_aile pour tout

# ==========================================
# 2. LE MODÈLE PHYSIQUE (FONCTION DÉRIVÉE)
# ==========================================
def modele_planeur(vecteur_etat, t):
    """
    Entrée : vecteur_etat = [x, y, vx, vy]
    Sortie : [vx, vy, ax, ay]
    """
    
    # A. Lecture de l'état actuel
    # On ne se sert pas de x et y pour le calcul des forces (sauf si l'atmosphère changeait avec l'altitude)
    vx = vecteur_etat[2]
    vy = vecteur_etat[3]
    
    # B. Calcul des vitesses et forces
    v_carre = vx**2 + vy**2
    v_norme = np.sqrt(v_carre)
    
    # Forces Aérodynamiques (Normes)
    force_portance = k_portance * v_carre
    force_trainee  = k_trainee  * v_carre
    
    # C. Projection des forces (Le moment délicat !)
    # Le vecteur vitesse indique la direction du mouvement.
    # La traînée est OPPOSÉE à la vitesse.
    # La portance est PERPENDICULAIRE à la vitesse (vers le "haut" du planeur).
    
    # Cosinus et Sinus de l'angle de la trajectoire
    # (Astuce : cos = vx/v, sin = vy/v pour éviter les calculs d'angles)
    cos_alpha = vx / v_norme
    sin_alpha = vy / v_norme
    
    # Application du PFD (Somme des forces = m * a) sur les axes X et Y
    
    # Sur l'axe X (Horizontal) :
    # - La traînée freine (s'oppose à vx) -> -T * cos
    # - La portance tire vers l'arrière quand on monte/descend -> -P * sin
    fx = -force_trainee * cos_alpha - force_portance * sin_alpha
    ax = fx / m
    
    # Sur l'axe Y (Vertical) :
    # - La traînée freine (s'oppose à vy) -> -T * sin
    # - La portance soutient l'avion -> +P * cos
    # - Le poids tire vers le bas -> -m*g
    fy = -force_trainee * sin_alpha + force_portance * cos_alpha - m*g
    ay = fy / m
    
    # D. Retour des dérivées : [vitesse_x, vitesse_y, accel_x, accel_y]
    return [vx, vy, ax, ay]

# ==========================================
# 3. PRÉPARATION DE LA SIMULATION
# ==========================================

# A. Conditions Initiales
x0 = 0.0       # Position départ
y0 = 1000.0    # Altitude départ (m)
vx0 = 25.0     # Vitesse horizontale (m/s)
vy0 = 0.0      # Vitesse verticale (chute initiale nulle)

etat_initial = [x0, y0, vx0, vy0]

# B. Temps de simulation
t = np.linspace(0, 120, 1000) # 2 minutes de vol

# ==========================================
# 4. RÉSOLUTION NUMÉRIQUE
# ==========================================
solution = odeint(modele_planeur, etat_initial, t)

# Extraction des résultats pour la lisibilité
x_traj = solution[:, 0]
y_traj = solution[:, 1]
# vx_traj = solution[:, 2] # Pas utilisé pour le tracé principal
# vy_traj = solution[:, 3]

# ==========================================
# 5. VISUALISATION SCIENTIFIQUE
# ==========================================
plt.figure(figsize=(10, 6))

plt.plot(x_traj, y_traj, label="Trajectoire du planeur", color='blue')

plt.title(f"Simulation de vol plané (Départ à {y0}m)")
plt.xlabel("Distance parcourue au sol (m)")
plt.ylabel("Altitude (m)")
plt.axhline(0, color='black', linewidth=2) # Ligne du sol
plt.grid(True)
plt.legend()

# On force un ratio égal pour que la pente soit réaliste visuellement
plt.axis('equal') 

plt.show()