SciPy : Calcul Scientifique Avancé en Python

SciPy est la bibliothèque de calcul scientifique qui prolonge NumPy avec des algorithmes numériques éprouvés : optimisation, ajustement de courbes, statistiques, intégration d’équations différentielles, interpolation, traitement du signal et matrices creuses. Là où NumPy fournit la structure de données (le ndarray) et les opérations élémentaires, SciPy fournit les méthodes numériques de haut niveau, pour la plupart adossées à des bibliothèques Fortran et C de référence (LAPACK, MINPACK, FFTPACK). Ce guide parcourt les sous-modules les plus utilisés avec des exemples complets et exécutables.

SciPy par rapport à NumPy

La frontière est simple à retenir :

• NumPy : le tableau ndarray, l’arithmétique vectorisée, l’algèbre linéaire de base, la génération aléatoire.
• SciPy : les algorithmes qui opèrent sur ces tableaux, avec des implémentations robustes issues de décennies de calcul scientifique.

SciPy dépend de NumPy et s’importe sous-module par sous-module :

import numpy as np
from scipy import optimize, stats, integrate, interpolate, signal, sparse

# Vérifier la version installée
import scipy
print(scipy.__version__)     # 1.14.1 (par exemple)

Quelques recouvrements existent (scipy.linalg étend numpy.linalg avec des décompositions supplémentaires ; scipy.fft remplace avantageusement numpy.fft), et dans ces cas la version SciPy est généralement plus complète.

scipy.optimize : minimisation et ajustement de courbes

Minimiser une fonction avec minimize

minimize trouve un minimum local d’une fonction scalaire de plusieurs variables.

from scipy.optimize import minimize

# Fonction de Rosenbrock, un banc d'essai classique
def rosenbrock(p):
    x, y = p
    return (1 - x) ** 2 + 100 * (y - x ** 2) ** 2

resultat = minimize(rosenbrock, x0=[0.0, 0.0], method="Nelder-Mead")

print(resultat.success)      # True
print(resultat.x)            # [0.99998736 0.99997369] -> minimum en (1, 1)
print(resultat.fun)          # 5.9e-10 (valeur au minimum, quasi nulle)
print(resultat.nit)          # nombre d'itérations

Le choix de la méthode dépend du problème : "Nelder-Mead" ne requiert pas de gradient, "BFGS" (défaut sans contraintes) l’estime par différences finies, "L-BFGS-B" accepte des bornes, "SLSQP" gère des contraintes générales.

# Minimisation avec bornes : x dans [0, 2], y dans [0, 2]
res = minimize(rosenbrock, x0=[0.5, 0.5], method="L-BFGS-B",
               bounds=[(0, 2), (0, 2)])
print(res.x)                 # [1. 1.]

Pour trouver la racine d’une équation plutôt qu’un minimum, utilisez optimize.root_scalar ou optimize.brentq.

Ajuster un modèle à des données avec curve_fit

curve_fit estime les paramètres d’un modèle par moindres carrés non linéaires. Exemple : ajuster une décroissance exponentielle à des mesures bruitées.

from scipy.optimize import curve_fit

# Modèle : la première variable est toujours x, les suivantes sont les paramètres
def modele(x, a, tau, c):
    return a * np.exp(-x / tau) + c

# Données synthétiques bruitées
rng = np.random.default_rng(42)
x_donnees = np.linspace(0, 10, 50)
y_vrai = modele(x_donnees, a=3.0, tau=2.5, c=0.5)
y_donnees = y_vrai + rng.normal(0, 0.15, size=x_donnees.size)

# Ajustement (p0 : estimation initiale des paramètres)
params, covariance = curve_fit(modele, x_donnees, y_donnees, p0=[1, 1, 0])

a_est, tau_est, c_est = params
incertitudes = np.sqrt(np.diag(covariance))

print(f"a   = {a_est:.3f} ± {incertitudes[0]:.3f}")   # a   = 3.011 ± 0.052
print(f"tau = {tau_est:.3f} ± {incertitudes[1]:.3f}") # tau = 2.487 ± 0.098
print(f"c   = {c_est:.3f} ± {incertitudes[2]:.3f}")   # c   = 0.503 ± 0.031

La diagonale de la matrice de covariance donne la variance de chaque paramètre estimé : sa racine carrée fournit une incertitude à un écart-type. Fournissez toujours un p0 raisonnable pour les modèles non linéaires, sous peine de convergence vers un mauvais optimum.

scipy.stats : distributions et tests statistiques

Travailler avec des distributions

Chaque distribution expose la même interface : pdf (densité), cdf (fonction de répartition), ppf (quantiles, inverse de la cdf) et rvs (échantillonnage).

from scipy import stats

# Loi normale de moyenne 100 et d'écart-type 15
loi = stats.norm(loc=100, scale=15)

print(loi.pdf(100))          # 0.02659... (densité au sommet)
print(loi.cdf(115))          # 0.8413... (P(X <= 115))
print(loi.ppf(0.975))        # 129.399... (quantile à 97,5 %)
print(loi.rvs(size=3, random_state=42))   # [107.45  84.44 109.72]

# Autres lois du même moule : stats.t, stats.chi2, stats.binom,
# stats.poisson, stats.expon, stats.uniform...

Test t : comparer deux moyennes

rng = np.random.default_rng(1)
groupe_a = rng.normal(loc=52, scale=8, size=40)   # ex. scores avec méthode A
groupe_b = rng.normal(loc=48, scale=8, size=40)   # ex. scores avec méthode B

t_stat, p_valeur = stats.ttest_ind(groupe_a, groupe_b)
print(f"t = {t_stat:.3f}, p = {p_valeur:.4f}")
# t = 2.279, p = 0.0255 (par exemple)

if p_valeur < 0.05:
    print("Différence significative au seuil de 5 %")
# Variante appariée (mesures avant/après sur les mêmes sujets) : stats.ttest_rel
# Variances inégales : ttest_ind(..., equal_var=False) (test de Welch)

Test du khi-deux : indépendance de deux variables catégorielles

# Tableau de contingence : lignes = formation suivie, colonnes = réussite
observations = np.array([[45, 15],    # avec formation : 45 réussites, 15 échecs
                         [30, 30]])   # sans formation : 30 réussites, 30 échecs

chi2, p, ddl, attendus = stats.chi2_contingency(observations)
print(f"chi2 = {chi2:.3f}, p = {p:.4f}, ddl = {ddl}")
# chi2 = 6.806, p = 0.0091, ddl = 1
print(attendus)              # effectifs attendus sous l'hypothèse d'indépendance
# [[37.5 22.5]
#  [37.5 22.5]]

Corrélations

x = rng.normal(size=100)
y = 0.7 * x + rng.normal(scale=0.5, size=100)

r, p = stats.pearsonr(x, y)          # corrélation linéaire
print(f"Pearson  r = {r:.3f}, p = {p:.2e}")   # r = 0.812, p = 1.5e-24

rho, p = stats.spearmanr(x, y)       # corrélation de rangs (relations monotones)
print(f"Spearman rho = {rho:.3f}")   # rho = 0.795

Utilisez Spearman quand la relation est monotone mais non linéaire, ou en présence de valeurs aberrantes.

scipy.integrate : intégrales et équations différentielles

Intégration numérique avec quad

from scipy.integrate import quad

# Intégrale de exp(-x²) entre 0 et l'infini (vaut sqrt(pi)/2)
valeur, erreur = quad(lambda x: np.exp(-x ** 2), 0, np.inf)
print(valeur)                        # 0.8862269254527579
print(np.sqrt(np.pi) / 2)            # 0.8862269254527579
print(erreur)                        # 7.1e-09 (estimation de l'erreur)

Résoudre une EDO avec solve_ivp

Exemple : l’oscillateur harmonique amorti, décrit par l’équation du second ordre x” + 2ζω x’ + ω² x = 0, réécrite en système du premier ordre.

from scipy.integrate import solve_ivp

omega = 2.0    # pulsation propre
zeta = 0.15    # taux d'amortissement

def oscillateur(t, y):
    x, v = y                              # position, vitesse
    return [v, -2 * zeta * omega * v - omega ** 2 * x]

solution = solve_ivp(
    oscillateur,
    t_span=(0, 20),                       # intervalle d'intégration
    y0=[1.0, 0.0],                        # position 1, vitesse nulle à t=0
    t_eval=np.linspace(0, 20, 400),       # points où évaluer la solution
    method="RK45",                        # Runge-Kutta adaptatif (défaut)
)

print(solution.success)                   # True
print(solution.y.shape)                   # (2, 400) -> position et vitesse
print(solution.y[0, :3])                  # [1.  0.995 0.980] (position amortie)

# Tracé de la position au cours du temps
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(solution.t, solution.y[0])
ax.set_xlabel("Temps (s)")
ax.set_ylabel("Position")
ax.set_title("Oscillateur harmonique amorti")
plt.show()

Pour les systèmes raides (constantes de temps très différentes), passez method="Radau" ou method="BDF".

scipy.interpolate : interpolation de données

L’interpolation construit une fonction continue passant par des points de mesure discrets.

from scipy.interpolate import CubicSpline, interp1d

# Points de mesure épars
x_mesures = np.array([0, 1, 2, 3, 4, 5], dtype=float)
y_mesures = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])

# Spline cubique : lisse, dérivable, recommandée
spline = CubicSpline(x_mesures, y_mesures)

x_fin = np.linspace(0, 5, 100)
y_interpole = spline(x_fin)

print(spline(2.5))                   # 0.5546... (valeur interpolée entre 2 et 3)
print(spline(2.5, 1))                # dérivée première au même point

# Interpolation linéaire simple
lineaire = interp1d(x_mesures, y_mesures, kind="linear")
print(lineaire(2.5))                 # 0.5 (moyenne des points voisins)

Attention : l’interpolation n’est fiable qu’à l’intérieur du domaine des mesures. Pour extrapoler ou pour lisser des données bruitées, préférez un ajustement de modèle (curve_fit) ou une spline de lissage (scipy.interpolate.make_smoothing_spline).

scipy.signal et scipy.fft : traitement du signal

Analyse spectrale avec la FFT

from scipy.fft import fft, fftfreq

# Signal échantillonné : somme de deux sinusoïdes (5 Hz et 50 Hz) plus du bruit
fe = 500                              # fréquence d'échantillonnage (Hz)
duree = 2.0
t = np.linspace(0, duree, int(fe * duree), endpoint=False)
rng = np.random.default_rng(7)
signal_brut = (np.sin(2 * np.pi * 5 * t)
               + 0.5 * np.sin(2 * np.pi * 50 * t)
               + 0.3 * rng.normal(size=t.size))

# Transformée de Fourier
spectre = fft(signal_brut)
frequences = fftfreq(t.size, d=1 / fe)

# Amplitudes sur les fréquences positives
moitie = t.size // 2
amplitudes = 2 / t.size * np.abs(spectre[:moitie])

pics = frequences[:moitie][amplitudes > 0.3]
print(pics)                           # [ 5. 50.] -> les deux fréquences retrouvées

Filtrage avec un filtre de Butterworth

from scipy.signal import butter, filtfilt

# Filtre passe-bas d'ordre 4, fréquence de coupure 20 Hz
b, a = butter(N=4, Wn=20, btype="low", fs=fe)

# filtfilt applique le filtre dans les deux sens : aucun déphasage
signal_filtre = filtfilt(b, a, signal_brut)

# La composante à 50 Hz et le bruit haute fréquence sont supprimés,
# la sinusoïde à 5 Hz est préservée.

scipy.signal fournit aussi la détection de pics (find_peaks), les spectrogrammes (ShortTimeFFT), la convolution et la conception de nombreux types de filtres (Chebyshev, elliptiques, FIR).

from scipy.signal import find_peaks

pics_indices, proprietes = find_peaks(signal_filtre, height=0.5, distance=fe // 10)
print(len(pics_indices))              # 10 -> un pic par période de la sinusoïde 5 Hz

scipy.sparse : matrices creuses

Quand une matrice contient une écrasante majorité de zéros (graphes, systèmes d’équations issus de maillages, matrices termes-documents), la stocker densément gaspille mémoire et calcul.

from scipy import sparse

# Matrice creuse au format COO (construction) puis CSR (calcul)
lignes = np.array([0, 1, 2, 2])
colonnes = np.array([1, 2, 0, 2])
valeurs = np.array([3.0, 4.0, 5.0, 6.0])

m = sparse.coo_array((valeurs, (lignes, colonnes)), shape=(3, 3)).tocsr()
print(m.toarray())
# [[0. 3. 0.]
#  [0. 0. 4.]
#  [5. 0. 6.]]

print(m.nnz)                          # 4 éléments non nuls stockés (sur 9)

# Produit matrice-vecteur : seuls les éléments non nuls sont parcourus
v = np.array([1.0, 2.0, 3.0])
print(m @ v)                          # [ 6. 12. 23.]

# Résolution d'un système linéaire creux
from scipy.sparse.linalg import spsolve
b = np.array([6.0, 12.0, 23.0])
x = spsolve(m.tocsc(), b)
print(x)                              # [1. 2. 3.]

Ordres de grandeur : une matrice dense de 100 000 × 100 000 en float64 occuperait 80 Go ; la même matrice avec 0,01 % d’éléments non nuls tient en quelques mégaoctets au format CSR.

Tableau récapitulatif des sous-modules

Sous-module	Domaine	Fonctions phares
scipy.optimize	Optimisation, ajustement, racines	minimize, curve_fit, root_scalar, linprog
scipy.stats	Distributions, tests, corrélations	norm, ttest_ind, chi2_contingency, pearsonr
scipy.integrate	Intégrales, équations différentielles	quad, solve_ivp
scipy.interpolate	Interpolation, splines	CubicSpline, interp1d, make_smoothing_spline
scipy.signal	Filtrage, pics, spectrogrammes	butter, filtfilt, find_peaks
scipy.fft	Transformées de Fourier	fft, ifft, fftfreq, rfft
scipy.sparse	Matrices creuses	csr_array, coo_array, sparse.linalg.spsolve
scipy.linalg	Algèbre linéaire étendue	lu, svd, expm, solve_banded
scipy.spatial	Géométrie, distances	KDTree, distance.cdist, ConvexHull
scipy.ndimage	Traitement d’images n-D	gaussian_filter, label, rotate

Bonnes pratiques

• Importez les sous-modules explicitement (from scipy import optimize) : import scipy seul ne charge pas les sous-modules.
• Vérifiez toujours resultat.success après minimize ou solve_ivp : un algorithme qui n’a pas convergé renvoie quand même un résultat.
• Fournissez une estimation initiale p0 sensée à curve_fit : les moindres carrés non linéaires convergent vers un optimum local dépendant du point de départ.
• Interprétez les p-valeurs avec prudence : une p-valeur sous 0,05 ne mesure ni la taille de l’effet ni son importance pratique ; rapportez aussi les moyennes et intervalles de confiance.
• Choisissez la méthode adaptée au problème : Radau ou BDF pour les EDO raides, L-BFGS-B dès qu’il y a des bornes, le test de Welch quand les variances diffèrent.
• N’extrapolez pas avec une interpolation : hors du domaine des mesures, une spline cubique diverge rapidement.
• Passez aux matrices creuses dès que la densité d’éléments non nuls descend sous quelques pour cent sur de grandes matrices.

Conclusion

SciPy transforme Python en environnement de calcul scientifique complet : optimisation et ajustement de modèles avec optimize, inférence statistique avec stats, résolution d’équations différentielles avec integrate, reconstruction de fonctions avec interpolate, analyse fréquentielle et filtrage avec signal et fft, et passage à l’échelle avec sparse. Chaque sous-module encapsule des algorithmes de référence validés depuis des décennies, ce qui vous évite de réimplémenter des méthodes numériques délicates. Combiné à NumPy pour les structures de données, Pandas pour les tableaux étiquetés et Matplotlib pour la visualisation, SciPy complète la boîte à outils standard de tout travail scientifique ou d’ingénierie en Python.