Résoudre une Equation Différentielle Ordinaire (EDO) en Python avec SciPy — exemple sur un capteur de givre aéronautique
Dans ma carrière d’ingénieur aéronautique, je n’ai jamais résolu une équation différentielle à la main sur un vrai projet. Les logiciels métier — SAMCEF, ABAQUS, NASTRAN — s’en occupaient. On définissait la géométrie, les matériaux, les conditions aux limites et le logiciel calculait les fréquences propres et les modes.
Mais comprendre ce que le logiciel fait sous le capot change complètement la façon dont on interprète les résultats. Et pouvoir reproduire un calcul de dynamique simple en Python, sans logiciel métier, ouvre des possibilités que les grands outils ne permettent pas ; explorer rapidement des variantes de design, visualiser l’effet d’un paramètre, automatiser une étude paramétrique.
Mon projet de fin d’études m’a confronté exactement à ce type de problème. Il portait sur un capteur de givre aéronautique, un dispositif dont le principe de fonctionnement repose entièrement sur la dynamique vibratoire. Quand du givre se dépose sur le capteur, sa masse augmente et sa fréquence propre varie. En mesurant ce décalage de fréquence, on détecte le givrage en temps réel.
C’est ce problème qu’on va résoudre aujourd’hui avec Python et SciPy.
Le principe physique du capteur de givre
Un capteur de givre vibrant peut être modélisé, en première approximation, comme un système masse-ressort-amortisseur — le modèle de base de la dynamique des structures.
m·ẍ + c·ẋ + k·x = F(t)
Avec :
m : masse totale (capteur + givre déposé) [kg]
c : coefficient d'amortissement [N·s/m]
k : raideur du capteur [N/m]
x : déplacement [m]
F : force d'excitation [N]La fréquence propre non amortie du système est :
f₀ = (1/2π) × √(k/m)Ce que ça dit physiquement : si la raideur k est fixe (le capteur ne change pas) et que la masse m augmente (du givre se dépose), alors f₀ diminue. C’est ce signal de décalage de fréquence que le système de détection mesure.
L’objectif du projet de fin d’études : simuler le comportement dynamique du capteur pour différentes quantités de givre et identifier les changements de design qui maximisent la sensibilité de détection.
Transformer une EDO d’ordre 2 en système d’ordre 1
SciPy résout des systèmes d’EDO du premier ordre. Notre équation est du second ordre. On la transforme en introduisant deux variables d’état :
y₁ = x (déplacement)
y₂ = ẋ (vitesse)
Alors :
ẏ₁ = y₂
ẏ₂ = (F(t) - c·y₂ - k·y₁) / mC’est la forme que SciPy attend — un vecteur de dérivées à chaque instant.
Le code Python complet
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from numpy.fft import fft, fftfreq
# ============================================================
# PARAMÈTRES DU CAPTEUR DE GIVRE
# ============================================================
K = 1500.0 # Raideur [N/m]
C = 0.8 # Amortissement [N·s/m]
M_CAPTEUR = 0.005 # Masse du capteur seul [kg]
# Scénarios de givrage simulés
SCENARIOS_GIVRE = {
'Sans givre' : 0.000,
'Givrage léger' : 0.001,
'Givrage modéré': 0.003,
'Givrage fort' : 0.006,
}
# Conditions initiales
X0 = 0.001 # déplacement initial [m]
V0 = 0.0 # vitesse initiale [m/s]
# Durée de simulation
T_FIN = 0.5
T_EVAL = np.linspace(0, T_FIN, 5000)
# ============================================================
# DÉFINITION DU SYSTÈME D'EDO
# ============================================================
def systeme_capteur(t, y, m, c, k):
"""
Système masse-ressort-amortisseur — réponse libre.
y[0] = x (déplacement)
y[1] = ẋ (vitesse)
Retourne [ẋ, ẍ] — les dérivées des deux variables d'état.
"""
x, v = y
dxdt = v
dvdt = (-c * v - k * x) / m
return [dxdt, dvdt]
# ============================================================
# RÉSOLUTION POUR CHAQUE SCÉNARIO
# ============================================================
resultats = {}
for scenario, masse_givre in SCENARIOS_GIVRE.items():
m_total = M_CAPTEUR + masse_givre
f0 = (1 / (2 * np.pi)) * np.sqrt(K / m_total)
solution = solve_ivp(
fun=systeme_capteur,
t_span=[0, T_FIN],
y0=[X0, V0],
t_eval=T_EVAL,
args=(m_total, C, K),
method='RK45',
rtol=1e-8,
atol=1e-10,
)
resultats[scenario] = {
'temps' : solution.t,
'deplacement' : solution.y[0],
'vitesse' : solution.y[1],
'm_total' : m_total,
'f0' : f0,
'masse_givre' : masse_givre,
}
f0_ref = (1 / (2 * np.pi)) * np.sqrt(K / M_CAPTEUR)
print(f"{scenario:20s} | m = {m_total*1000:.1f} g "
f"| f₀ = {f0:.2f} Hz "
f"| Δf = {f0 - f0_ref:.2f} Hz")
# ============================================================
# ANALYSE FRÉQUENTIELLE PAR FFT
# ============================================================
print("nAnalyse fréquentielle par FFT :")
for scenario, r in resultats.items():
signal = r['deplacement']
dt = T_EVAL[1] - T_EVAL[0]
n = len(signal)
spectre = np.abs(fft(signal))[:n//2]
freqs = fftfreq(n, dt)[:n//2]
f_mesuree = freqs[np.argmax(spectre)]
print(f"{scenario:20s} | f₀ théorique = {r['f0']:.2f} Hz "
f"| f₀ FFT = {f_mesuree:.2f} Hz")
# ============================================================
# VISUALISATION
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle(
'Capteur de givre — Analyse dynamique par scénario de givrage',
fontsize=13, fontweight='bold'
)
couleurs = ['#27ae60', '#f5a623', '#e94560', '#8e44ad']
# Graphique 1 : Réponses temporelles
ax1 = axes[0, 0]
for (scenario, r), couleur in zip(resultats.items(), couleurs):
ax1.plot(r['temps'] * 1000, r['deplacement'] * 1000,
color=couleur, lw=1.5,
label=f"{scenario} (f₀={r['f0']:.1f} Hz)",
alpha=0.85)
ax1.set_xlabel('Temps (ms)')
ax1.set_ylabel('Déplacement (mm)')
ax1.set_title('Réponse libre — comparaison des scénarios')
ax1.legend(fontsize=8)
ax1.grid(True, alpha=0.3)
# Graphique 2 : Zoom 50 ms
ax2 = axes[0, 1]
idx_zoom = T_EVAL < 0.05
for (scenario, r), couleur in zip(resultats.items(), couleurs):
ax2.plot(r['temps'][idx_zoom] * 1000,
r['deplacement'][idx_zoom] * 1000,
color=couleur, lw=1.8, label=scenario, alpha=0.85)
ax2.set_xlabel('Temps (ms)')
ax2.set_ylabel('Déplacement (mm)')
ax2.set_title('Zoom — premiers 50 msn(décalage de phase visible)')
ax2.legend(fontsize=8)
ax2.grid(True, alpha=0.3)
# Graphique 3 : Fréquence propre vs masse de givre
ax3 = axes[1, 0]
masses = [r['masse_givre'] * 1000 for r in resultats.values()]
f0s = [r['f0'] for r in resultats.values()]
ax3.plot(masses, f0s, 'o-', color='#e94560', lw=2, ms=8, zorder=5)
for m, f in zip(masses, f0s):
ax3.annotate(f'{f:.1f} Hz',
xy=(m, f),
xytext=(m + 0.1, f + 0.3),
fontsize=8, color='#333333')
ax3.set_xlabel('Masse de givre (g)')
ax3.set_ylabel('Fréquence propre f₀ (Hz)')
ax3.set_title('Sensibilité du capteurn(décalage f₀ en fonction du givrage)')
ax3.grid(True, alpha=0.3)
# Graphique 4 : Spectre FFT
ax4 = axes[1, 1]
for (scenario, r), couleur in zip(resultats.items(), couleurs):
signal = r['deplacement']
dt = T_EVAL[1] - T_EVAL[0]
n = len(signal)
spectre = np.abs(fft(signal))[:n//2]
freqs = fftfreq(n, dt)[:n//2]
spectre_norm = spectre / spectre.max()
idx_freq = freqs < 200
ax4.plot(freqs[idx_freq], spectre_norm[idx_freq],
color=couleur, lw=1.5, label=scenario, alpha=0.85)
ax4.set_xlabel('Fréquence (Hz)')
ax4.set_ylabel('Amplitude normalisée')
ax4.set_title('Spectre FFT — pic de résonance par scénarion'
'(déplacement du pic = détection du givre)')
ax4.legend(fontsize=8)
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('capteur_givre_dynamique.png', dpi=150, bbox_inches='tight')
plt.show()
Ce que les graphiques nous disent
Graphique 1 — Réponses temporelles : les quatre courbes ont la même enveloppe d’amortissement (même C et même K) mais des périodes différentes. Plus il y a de givre, plus la période est longue — la fréquence d’oscillation diminue.
Graphique 2 — Zoom sur 50 ms : le décalage de phase entre les courbes devient visible dès les premiers cycles. C’est ce décalage que le système électronique du capteur mesure en temps réel.
Graphique 3 — Sensibilité : c’est le graphique clé pour le design. Il montre comment la fréquence propre varie en fonction de la masse de givre. La pente de cette courbe est la sensibilité du capteur — plus elle est prononcée, plus le capteur est précis.
Graphique 4 — Spectre FFT (Transformée de Fourrier) : les pics de résonance se déplacent vers la gauche avec l’accumulation de givre. C’est la signature fréquentielle du givrage — celle que le système de détection lit en vol.
L’étude paramétrique — l’intérêt réel de Python
Ce que Python permet et qu’un logiciel métier rend fastidieux : lancer automatiquement des dizaines de variantes de design et comparer les résultats.
# Quelle raideur maximise la sensibilité de détection ?
raideurs = np.linspace(500, 5000, 20) # N/m
masse_givre_ref = 0.003 # 3g de givre
sensibilites = []
for k in raideurs:
f0_sans = (1/(2*np.pi)) * np.sqrt(k / M_CAPTEUR)
f0_avec = (1/(2*np.pi)) * np.sqrt(k / (M_CAPTEUR + masse_givre_ref))
delta_f = f0_sans - f0_avec
sensibilites.append(delta_f)
k_optimal = raideurs[np.argmax(sensibilites)]
print(f"Raideur optimale : {k_optimal:.0f} N/m")
print(f"Décalage de fréquence maximal : {max(sensibilites):.2f} Hz")
plt.figure(figsize=(8, 5))
plt.plot(raideurs, sensibilites, color='#e94560', lw=2)
plt.axvline(k_optimal, color='#f5a623', ls='--', lw=1.5,
label=f'k optimal = {k_optimal:.0f} N/m')
plt.xlabel('Raideur k (N/m)')
plt.ylabel('Décalage de fréquence Δf (Hz)')
plt.title('Optimisation du capteur de givren'
'Sensibilité en fonction de la raideur')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('optimisation_capteur.png', dpi=150)
plt.show()
C’est ce type d’étude paramétrique que mon projet de fin d’études réalisait — trouver le design qui maximise la sensibilité de détection pour une masse de givre donnée. Avec un logiciel métier, chaque variante demandait de reconfigurer un modèle. Avec Python, on boucle sur 20 variantes en quelques secondes.
« Dans la réalité, la raideur d’un capteur de givre industriel
ne se calcule pas par une formule analytique simple — la géométrie
est bien plus complexe qu’une poutre encastrée. C’est précisément
pour ça qu’on utilise des logiciels éléments finis comme Samcef
pour l’identifier. Python intervient ensuite : une fois k et m
connus depuis le modèle EF, on peut simuler le comportement
dynamique, automatiser l’étude paramétrique, et explorer
l’espace de design rapidement. »
Ce que SciPy fait sous le capot
La méthode RK45 utilisée par solve_ivp est un solveur de Runge-Kutta d’ordre 4-5 adaptatif. Concrètement :
Il avance dans le temps pas à pas, en estimant la dérivée à plusieurs points intermédiaires pour construire une solution précise. Le pas de temps est ajusté automatiquement — petit quand le signal varie vite (début de la réponse), plus grand quand il varie lentement (fin de l’amortissement). Les paramètres rtol et atol contrôlent la précision : plus ils sont petits, plus la solution est précise mais plus le calcul est long.
C’est ce que votre logiciel métier (Samcef, Abaqus) fait aussi — avec des méthodes numériques bien plus sophistiquées pour des systèmes bien plus complexes. Python ne remplace pas ces outils pour des structures réelles. Mais il permet de comprendre, d’explorer, et d’automatiser des analyses rapides sur des modèles simplifiés.
Conclusion
Résoudre une EDO en Python ne demande pas d’être mathématicien. Il faut comprendre deux choses : comment transformer son équation en système du premier ordre, et comment appeler solve_ivp avec les bons arguments.
Le reste — le choix du solveur, le contrôle du pas de temps, la gestion de la précision — SciPy s’en occupe.
Et une fois que la résolution est automatisée, l’étude paramétrique devient triviale. C’est là que Python apporte sa vraie valeur pour un ingénieur mécanicien : pas dans le calcul d’un cas unique, mais dans l’exploration rapide de l’espace de design.
🚀 Vous souhaitez aller plus loin ?
Téléchargez le Pack Complet « De Facto »
J’ai condensé une méthode dans le ‘De Facto Data Automation Kit’:
• 📊 des données capteurs brutes utilisées dans le challenge 7 jours ;
• 🐍 le script Python de nettoyage ;
• ✅ la checklist de contrôle des données ;
• 📈 l’exemple avant/après ;
• 🔧 les scripts des 7 étapes du challenge ;
• 📄 le guide complet.
En attendant, j’aimerais connaître votre réalité de terrain
💬 Vous avez utilisé des méthodes numériques dans votre parcours d’ingénieur ? À la main, avec MATLAB, ou avec autre chose ? Dites-moi en commentaire.
Laisser un commentaire