Un lien vers Google Colab qui peut faire tourner le notebook qui contient les exercices se trouve ici :
https://colab.research.google.com/github/nickcollins-craft/USMA1Q-Methodes-Numeriques/blob/main/
Plan du cours¶
| Jour | Sujet |
|---|---|
| 03 mars | Fondamentaux de la programmation en Python |
| 04 mars (matin) | Fonctions, listes approfondies, fichiers et graphiques |
| 04 mars (après-midi) | Calcul numérique avec NumPy |
| 11 mars (matin) | Tableaux en mémoire, matrices & précision numérique |
| 11 mars (après-midi) | Résolution des systèmes linéaires |
| 18 mars | Interpolation & recherche des racines |
| 25 mars | Intégration numérique |
| 02 avril | Équations différentielles |
| 08 avril | Examen |
Plan de la séance (indicatif)¶
| Horaire | Sujet |
|---|---|
| 0:00 – 1:00 | Méthodes directes : np.linalg.solve(), factorisation LU |
| 1:00 – 1:50 | Conditionnement & quand ça tourne mal |
| 1:50 – 2:00 | ☕ Pause |
| 2:00 – 2:40 | Moindres carrés & ajustement de courbes |
| 2:40 – 3:40 | Mini-projet noté — Introduction |
Motivation — pourquoi résoudre des systèmes linéaires ?¶
Beaucoup de problèmes de science des matériaux se réduisent à résoudre $Ax = b$ :
| Domaine | $A$ | $x$ | $b$ |
|---|---|---|---|
| Conduction thermique | matrice de conductivité | températures inconnues | conditions aux limites |
| Treillis (éléments finis) | matrice de rigidité | déplacements | vecteur de force |
| Diffusion | matrice de diffusivité | concentrations | flux imposés |
| Ajustement de courbes | matrice de Vandermonde | coefficients | valeurs mesurées |
Résoudre ce problème efficacement et précisément est une compétence de base en calcul scientifique.
Dans les séances futures, nous rencontrerons $Ax = b$ sous diverses formes déguisées.
1 · Méthodes directes pour $Ax = b$¶
Élimination de Gauss — le concept¶
L'élimination de Gauss (Gaussian elimination) transforme $Ax = b$ en un système triangulaire supérieur $Ux = c$ par opérations sur les lignes, puis remonte (back-substitution) pour trouver $x$.
[2 1 3 | 1] → [2 1 3 | 1] → [2 1 3 | 1]
[4 3 5 | 3] → [0 1 -1 | 1] → [0 1 -1 | 1]
[6 2 8 | 4] → [0 -1 -1 | 1] → [0 0 -2 | 2]
Coût : $O(n^3)$ opérations — pour $n = 1000$, c'est $10^9$ multiplications.
En pratique : on n'implémente jamais l'élimination de Gauss à la main.
On utilise np.linalg.solve(), qui fait mieux.
Factorisation LU¶
La factorisation LU (avec permutation partiel) est un algorithm qui arrive au meme but, sans le coût élévé de l'élimination de Gauss.
Elle décompose $A$ en : $PA = LU$
- $P$ : matrice de permutation (permutation matrix) — réordonne les lignes pour la stabilité numérique (pivotage)
- $L$ : matrice triangulaire inférieure (lower) avec des 1 sur la diagonale
- $U$ : matrice triangulaire supérieure (upper)
Si on doit résoudre $Ax = b_1, Ax = b_2, \ldots, Ax = b_k$ (même $A$, différents seconds membres), on factorise $A$ une seule fois puis on résout $k$ systèmes triangulaires — chacun coûte seulement $O(n^2)$.
from scipy.linalg import lu
P, L, U = lu(A) # lu = LU factorisation
# Vérification : P @ A == L @ U
En mécanique des structures : on factorise la matrice de rigidité une fois, puis on applique différents chargements.
np.linalg.solve() — l'outil pratique¶
import numpy as np
A = np.array([[3, 1],
[1, 2]], dtype=float)
b = np.array([9, 8], dtype=float)
x = np.linalg.solve(A, b) # solve = résoudre
# → [2. 3.] car 3×2 + 1×3 = 9 et 1×2 + 2×3 = 8
Toujours vérifier la solution :
np.allclose(A @ x, b) # → True si la solution est correcte
Résolution pour plusieurs seconds membres :
B = np.array([[9, 6],
[8, 5]], dtype=float)
X = np.linalg.solve(A, B) # résout AX = B pour chaque colonne de B
⚠️ Ne jamais utiliser np.linalg.inv(A) @ b¶
# MAUVAIS — lent et moins précis numériquement
x = np.linalg.inv(A) @ b
# BON — utiliser solve()
x = np.linalg.solve(A, b)
Pourquoi ?
- Performance : calculer $A^{-1}$ coûte $O(n^3)$ et multiplier par $b$ coûte $O(n^2)$ supplémentaire.
solve()est plus court d'une étape. - Précision :
solve()utilise le pivotage partiel et évite les divisions par de très petits nombres. L'inversion explicite peut amplifier les erreurs d'arrondi. - Interprétabilité : dans les formules mathématiques, $A^{-1}b$ est courant — mais en code,
solve(A, b)exprime l'intention (résoudre un système) plus clairement.
# ── Démonstration : np.linalg.solve() et conduction thermique ─
import numpy as np
# Problème : barre d'acier, conduction thermique 1D stationnaire
# 5 nœuds intérieurs, T_gauche = 100 °C, T_droite = 25 °C
# Équation discrétisée : T[i - 1] - 2*T[i] + T[i + 1] = 0 (pas de source interne) → système tridiagonal
# (avec des entries non-nul seulement sur la diagonale principale et les diagonales adjacentes)
n = 5 # nombre de nœuds intérieurs
T_gauche = 100.0 # °C
T_droite = 25.0 # °C
# Matrice K (tridiagonale)
K = np.zeros((n, n))
for i in range(n):
K[i, i] = -2.0
if i > 0:
K[i, i - 1] = 1.0
if i < n-1:
K[i, i + 1] = 1.0
# Membres droits
b = np.zeros(n)
b[0] = -T_gauche # condition T_gauche = 100 °C
b[-1] = -T_droite # condition T_droite = 25 °C
print("Matrice K :", K)
print("Vecteur b :", b)
# Résolution
T_int = np.linalg.solve(K, b)
print("Températures intérieures (°C) :", np.round(T_int, 2))
# Vérification
print("K @ T_int == b :", np.allclose(K @ T_int, b))
# Profil complet (nœuds + frontières)
T_complet = np.concatenate([[T_gauche], T_int, [T_droite]])
x_noeuds = np.linspace(0, 1, n + 2)
print("Profil complet T(x) :", np.round(T_complet, 2))
À vous de faire l'exercice 1 !¶
2 · Conditionnement¶
Problème bien conditionné vs mal conditionné¶
Bien conditionné (well-conditioned) : une petite perturbation de $b$ produit une petite perturbation de $x$.
Mal conditionné (ill-conditioned) : une petite perturbation de $b$ produit une grande perturbation de $x$.
*Exemple géométrique :**
Système mal conditionné — droites presque parallèles :
$y = 2x + 1$ & $y = 2.001x + 0.999$
→ Les pentes diffèrent de 0.001 : un petit déplacement d'une droite change drastiquement l'intersection.
Système bien conditionné — droites à angle franc :
$y = 2x + 1$ & $y = 3x - 1$
→ Pentes 2 et 3 : l'intersection est stable face aux perturbations.
En termes pratiques : des mesures avec 1 % d'incertitude peuvent donner des résultats avec 1000 % d'incertitude si le système est mal conditionné.
Numéro de condition — np.linalg.cond(A)¶
Le numéro de condition (condition number) $\kappa(A)$ mesure la sensibilité de la solution aux perturbations :
$$\kappa(A) = \|A\| \cdot \|A^{-1}\|$$
Règle pratique :
| $\kappa(A)$ | Interprétation |
|---|---|
| ~1 | Parfaitement conditionné |
| ~10² | On perd ~2 chiffres significatifs |
| ~10⁸ | On perd ~8 chiffres sur 15 → douteux |
| ~10¹⁵ | On perd tous les chiffres → résultat sans sens |
| ≥ 10¹⁶ | Matrice numériquement singulière |
kappa = np.linalg.cond(A) # cond = condition number
# Si kappa ≈ 10^k, la solution perd environ k chiffres significatifs sur les 15-16 de float64
Résidu — comment vérifier sa solution¶
Même si np.linalg.solve() retourne un résultat sans erreur, vous pouvez toujours vérifier avec :
$$r = b - Ax$$
Si $\|r\|$ est grande, soit le système est mal conditionné, soit il y a un bug.
x = np.linalg.solve(A, b)
r = b - A @ x
print("Résidu ||r|| =", np.linalg.norm(r))
print("Résidu relatif ||r||/||b|| =", np.linalg.norm(r) / np.linalg.norm(b))
Bon résidu ≠ bonne solution si $\kappa$ est grand :
Le résidu peut être petit ($10^{-10}$) pendant que l'erreur sur $x$ est grande ($10^5$).
# Vérification complète
print("cond(A) =", np.linalg.cond(A))
print("||r|| =", np.linalg.norm(b - A @ x))
print("solution =", np.allclose(A @ x, b))
# Démonstration : conditionnement
# Matrice bien conditionnée (identité)
A_bon = np.array([[3, 1],
[1, 4]], dtype=float)
b = np.array([7, 9], dtype=float)
x_bon = np.linalg.solve(A_bon, b)
print("=== Système bien conditionné ===")
print("cond(A) = ", round(np.linalg.cond(A_bon), 2))
print("x = ", x_bon)
print("résidu ||r|| = ", np.linalg.norm(b - A_bon @ x_bon))
# Perturbation de b
b_pert = b + np.array([0.001, 0.0]) # 0.1 % de perturbation sur b[0]
x_pert_bon = np.linalg.solve(A_bon, b_pert)
print("Δx (perturb. 0.1%) = ", np.round(x_pert_bon - x_bon, 6))
# Matrice de Hilbert — célèbre pour son mauvais conditionnement
n = 6
H = np.array([[1/(i + j + 1) for j in range(n)] for i in range(n)])
b_H = np.ones(n)
x_H = np.linalg.solve(H, b_H)
print("=== Matrice de Hilbert ===")
print("cond(H) = ", round(np.linalg.cond(H), 2))
print("résidu ||r|| = ", round(np.linalg.norm(b_H - H @ x_H), 2))
b_H_pert = b_H + 0.001 * np.ones(n) # 0.1 % de perturbation
x_H_pert = np.linalg.solve(H, b_H_pert)
print("||Δx|| (perturb. 0.1%) = ", round(np.linalg.norm(x_H_pert - x_H), 2))
print("erreur amplifiée d'un facteur κ !")
À vous de faire l'exercice 2 !¶
☕ Pause — 10 minutes¶
3 · Moindres carrés & ajustement de courbes¶
Systèmes surdéterminés (overdetermined systems)¶
Un système surdéterminé a plus d'équations que d'inconnues — il n'y a généralement pas de solution exacte.
**Exemple : **ajuster un module de Young sur 50 points de traction :
- 50 équations : $\sigma_i = E \cdot \varepsilon_i + \sigma_0$
- 2 inconnues : $E$ et $\sigma_0$
- Système $50 \times 2$ — surdéterminé
Solution aux moindres carrés (least squares) : trouver le $x$ qui minimise $\|Ax - b\|^2$.
Équations normales (normal equations) : la solution satisfait $A^T A x = A^T b$.
En pratique : on n'utilise jamais les équations normales directement (mauvais conditionnement).
On utilisenp.linalg.lstsq()qui est plus stable numériquement.
np.linalg.lstsq() — usage pratique¶
x, residuals, rank, sv = np.linalg.lstsq(A, b, rcond=None)
# lstsq = least squares (moindres carrés)
# x : vecteur solution
# residuals : ||Ax - b||² pour chaque colonne (si A est full-rank)
# rank : rang de A
# sv : valeurs singulières de A (on peut utiliser pour diagnostiquer le conditionnement)
# rcond=None : utiliser la précision machine (recommandé)
Exemple — ajustement linéaire :
t = np.array([0, 1, 2, 3, 4], dtype=float)
y = np.array([2.1, 3.9, 6.2, 7.8, 10.1], dtype=float)
# Modèle : y = a*t + b → système A @ [a, b] = y.
# Les données "t" et "y" vont dans la matrice A et le vectuer b respectivement
A = np.column_stack([t, np.ones_like(t)]) # column_stack = empiler en colonnes
coeffs, _, _, _ = np.linalg.lstsq(A, y, rcond=None) # Nous donnes quelques variables (residuels, rang, valeurs singulières) le nom _ pour indiquer qu'on ne les utilisent pas, même si on est obligé de donner une variable sur la coté gauche pour prendre le valeur
a, b = coeffs
print("a = ", round(a, 3), " b = ", round(b, 3))
np.polyfit() pour les polynômes :
# Ajustement polynômial de degré deg sur les données (x, y)
coeffs = np.polyfit(x, y, deg=1) # retour les coefficients des puissances en list, de plus grande à la plus petite
y_ajuste = np.polyval(coeffs, x)
Loi puissance et linéarisation¶
Beaucoup de lois en science des matériaux sont de la forme :
$$y = A \cdot x^n$$
Cette loi non linéaire se linéarise par passage au logarithme :
$$\ln y = \ln A + n \ln x$$
En posant $Y = \ln y$, $X = \ln x$, $B = \ln A$, on obtient :
$$Y = n \cdot X + B$$
C'est un problème linéaire — parfait pour lstsq() ou polyfit().
Lois courantes en fluage :
- Loi de puissance en contrainte : $\dot{\varepsilon} = A \sigma^n$ — linéarisée en $\ln\dot\varepsilon = \ln A + n \ln\sigma$
- Loi puissance en temps : $\varepsilon = A t^n$ — linéarisée en $\ln\varepsilon = \ln A + n \ln t$
# Démonstration : ajustement loi puissance
import numpy as np
import matplotlib.pyplot as plt
# Données simulées : loi de fluage ε = A * t^n
# Valeurs vraies : A = 0.002, n = 0.35
np.random.seed(42)
t_data = np.array([10, 30, 100, 300, 1000, 3000, 10000], dtype=float) # heures
eps_vrai = 0.002 * t_data**0.35
eps_data = eps_vrai * (1 + 0.05 * np.random.randn(len(t_data))) # + 5 % bruit
# Linéarisation
ln_t = np.log(t_data)
ln_eps = np.log(eps_data)
# Matrice A pour le système linéaire : [n, ln(A)]
A_mat = np.column_stack([ln_t, np.ones_like(ln_t)])
coeffs, residus, rang, sv = np.linalg.lstsq(A_mat, ln_eps, rcond=None)
n_ajuste = coeffs[0]
A_ajuste = np.exp(coeffs[1])
print("Paramètres vrais : A = 0.002, n = 0.35")
print("Paramètres ajustés: A = %.4f, n = %.4f" % (A_ajuste, n_ajuste))
# Visualisation
t_courbe = np.linspace(t_data.min(), t_data.max(), 200)
eps_ajuste = A_ajuste * t_courbe**n_ajuste
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
# Échelle log-log — la relation linéarisée apparaît droite
axes[0].loglog(t_data, eps_data, 'o', label='données')
axes[0].loglog(t_courbe, eps_ajuste, '-', label='ajustement')
axes[0].set_xlabel('Temps (h)')
axes[0].set_ylabel('Déformation de fluage')
axes[0].set_title('Loi puissance — échelle log-log')
axes[0].legend()
axes[0].grid(True, which='both', alpha=0.3)
# Résidus
eps_predit = A_ajuste * t_data**n_ajuste
residus_plot = eps_data - eps_predit
axes[1].semilogx(t_data, residus_plot, 'o')
axes[1].axhline(0, color='k', ls='--')
axes[1].set_xlabel('Temps (h)')
axes[1].set_ylabel('Résidu (données − ajustement)')
axes[1].set_title('Analyse des résidus')
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
À vous de faire l'exercice 3 !¶
4 · Mini-projet noté — Analyse d'une série d'essais de traction¶
Contexte¶
Vous avez des données expérimentales pour 5 éprouvettes d'un acier inoxydable 316L.
Pour chaque éprouvette, vous disposez de :
- Un fichier CSV de force–déplacement (colonnes :
force_N,deplacement_mm) - Les dimensions : diamètre initial $d_0$ (mm) et longueur de jauge $L_0$ (mm)
Objectif : caractériser les propriétés mécaniques de cet alliage.
Données disponibles :
donnees_traction_acier_316L.csv: colonneseprouvette(1–5),force_N,deplacement_mmdimensions_eprouvettes_acier_316L.csv: colonneseprouvette,diametre_mm,longueur_mm
Tâche 1 — Chargement et calcul des courbes
- Charger les fichiers CSV avec
np.genfromtxt()ounp.loadtxt() - Calculer la contrainte d'ingénierie $\sigma = F / A_0$ et la déformation d'ingénierie $\varepsilon = \Delta L / L_0$
- Tracer les 5 courbes $\sigma(\varepsilon)$ sur un même graphique, avec légendes et axes annotés
Tâche 2 — Résistance à la traction et déformation à la rupture
- Pour chaque éprouvette, identifier la contrainte maximale (ultimate tensile strength, UTS) et la déformation à la rupture
- Calculer la moyenne, l'écart-type et le coefficient de variation (CV = écart-type/moyenne × 100 %)
Tâche 3 — Module de Young
- Dans le domaine élastique (déformation < 0,2 %), ajuster une droite $\sigma = E \varepsilon$ via
np.linalg.lstsq() - Calculer $E$ pour chaque éprouvette
- Comparer avec la valeur attendue pour 316L (~190–210 GPa)
Tâche 4 — Discussion
- Y a-t-il des éprouvettes anormales ? Justifiez
Tâche 5 — Propagation d'incertitude
- Estimez l'incertitude sur UTS en prenant $\delta F = 50$ N et $\delta d = 0{,}01$ mm
Règles de code :
- Dès qu'un tableu est établie, pas de boucles
forsur les éléments (utilisez les opérations vectorielles NumPy). Il est autorisé d'utiliser une boucle pour parcourir les éprouvettes. - Fonctions pour toute opération répétée (une fonction = un calcul)
- Axes annotés avec unités sur chaque graphique
- Variables nommées de façon descriptive
Format de rendu :
- Notebook Jupyter (
.ipynb) - Une fonction par calcul clé (
calculer_contrainte(),calculer_module_young(), etc.)
Rendu : samedi soir (avant minuit)
Barème indicatif :
| Tâche | Points |
|---|---|
| Tâche 1 : courbes et calculs | 5 |
| Tâche 2 : UTS et déformation rupture | 3 |
| Tâche 3 : module de Young | 5 |
| Tâche 4 : discussion | 2 |
| Tâche 5 : incertitudes | 3 |
| Clarté et style du code | 2 |
Points clés de la séance¶
| Concept | Ce qu'il faut retenir |
|---|---|
np.linalg.solve(A, b) |
La façon correcte de résoudre $Ax = b$ — toujours vérifier avec np.allclose(A@x, b) |
Jamais np.linalg.inv() |
Plus lent et moins précis que solve() |
| Factorisation LU | Réutilisable pour plusieurs seconds membres $b$ |
| Numéro de condition | np.linalg.cond(A) — si $> 10^8$, méfiance ; si $> 10^{15}$, résultat sans valeur |
| Résidu | r = b - A @ x — calculer np.linalg.norm(r) après résolution pour vérifier |
np.linalg.lstsq() |
Systèmes surdéterminés — minimise $\|Ax-b\|^2$ |
| Linéarisation | Loi puissance $y = Ax^n$ → $\ln y = \ln A + n\ln x$ → problème linéaire |
rcond=None |
Toujours passer rcond=None à lstsq() pour éviter l'avertissement |