Solving an Ordinary Differential Equation (ODE) in Python with SciPy — example on an aeronautical ice sensor
In my career as an aeronautical engineer, I never solved a differential equation by hand on a real project. The specialized software — SAMCEF, ABAQUS, NASTRAN — took care of that. We defined the geometry, the materials, the boundary conditions, and the software calculated the natural frequencies and modes.
But understanding what the software does under the hood completely changes how we interpret the results. And being able to reproduce a simple dynamics calculation in Python, without specialized software, opens up possibilities that large tools don’t allow: quickly exploring design variations, visualizing the effect of a parameter, automating a parametric study.
My final year project confronted me with precisely this type of problem. It focused on an aircraft icing sensor, a device whose operating principle relies entirely on vibrational dynamics. When icing forms on the sensor, its mass increases and its natural frequency changes. By measuring this frequency shift, icing can be detected in real time.
This is the problem we’re going to solve today with Python and SciPy.
The physical principle of the frost sensor
A vibrating frost sensor can be modeled, as a first approximation, as a mass-spring-damper system — the basic model of structural dynamics.
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]The undamped natural frequency of the system is:
f₀ = (1/2π) × √(k/m)What this means physically is: if the stiffness k is fixed (the sensor doesn’t change) and the mass m increases (frost is deposited), then f₀ decreases. It is this frequency shift signal that the detection system measures.
The objective of the final year project: to simulate the dynamic behavior of the sensor for different quantities of frost and to identify the design changes that maximize the detection sensitivity.
Transforming a second-order ODE into a first-order system
SciPy solves first-order ODE systems. Our equation is second-order. We transform it by introducing two state variables:
y₁ = x (déplacement)
y₂ = ẋ (vitesse)
Alors :
ẏ₁ = y₂
ẏ₂ = (F(t) - c·y₂ - k·y₁) / mThis is the form SciPy expects — a vector of derivatives at each instant.
The complete Python code
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()
What the graphs tell us
Figure 1 — Time responses: the four curves have the same damping envelope (same C and same K) but different periods. The more frost there is, the longer the period — the oscillation frequency decreases.
Figure 2 — Zooming in on 50 ms: the phase shift between the curves becomes visible from the very first cycles. This is the shift that the sensor’s electronic system measures in real time.
Graph 3 — Sensitivity: This is the key graph for the design. It shows how the natural frequency varies with the mass of frost. The slope of this curve represents the sensor’s sensitivity—the steeper the slope, the more accurate the sensor.
Figure 4 — FFT (Fourier Transform) Spectrum: the resonance peaks shift to the left with the accumulation of ice. This is the frequency signature of icing — the one that the detection system reads in flight.
Parametric analysis — the real value of Python
What Python allows, and what business software makes tedious: automatically launching dozens of design variants and comparing the results.
# 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()
This is the type of parametric study my final year project was conducting—finding the design that maximizes detection sensitivity for a given mass of frost. With specialized software, each variation required reconfiguring a model. With Python, you can loop through 20 variations in a few seconds.
“In reality, the stiffness of an industrial frost sensor
cannot be calculated using a simple analytical formula — the geometry
is far more complex than a cantilevered beam. This is precisely
why we use finite element software like Samcef
to identify it. Python then comes into play: once k and m
are known from the FE model, we can simulate the dynamic behavior, automate the parametric analysis, and quickly explore the design space.”
What SciPy does under the hood
The RK45 method used by solve_ivp is an adaptive 4th- or 5th-order Runge-Kutta solver. Specifically:
It progresses step by step through time, estimating the derivative at several intermediate points to construct an accurate solution. The time step is adjusted automatically—small when the signal changes rapidly (beginning of the response), larger when it changes slowly (end of damping). The parameters rtol and atol control the accuracy: the smaller they are, the more accurate the solution, but the longer the computation.
This is what your business software (Samcef, Abaqus) also does—with far more sophisticated numerical methods for far more complex systems. Python doesn’t replace these tools for real-world structures. But it allows you to understand, explore, and automate rapid analyses on simplified models.
Conclusion
Solving an ODE in Python doesn’t require being a mathematician. You need to understand two things: how to transform your equation into a first-order system, and how to call solve_ivp with the correct arguments.
The rest — the choice of solver, the control of the time step, the management of precision — SciPy takes care of.
And once the solution is automated, the parametric study becomes trivial. This is where Python brings its true value for a mechanical engineer: not in the calculation of a single case, but in the rapid exploration of the design space.
💬 Have you used numerical methods in your engineering career? By hand, with MATLAB, or with something else? Tell me in the comments.
📩 I am preparing a comprehensive training course on Python applied to engineering.
Leave a Reply