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 | Tableaux en mémoire : dtype, foulées, vues & copies |
| 1:00 – 1:50 | Construction et opérations matricielles |
| 1:50 – 2:00 | ☕ Pause |
| 2:00 – 2:40 | Précision numérique & virgule flottante |
| 2:40 – 3:45 | Analyse d'erreur & propagation d'incertitude |
1 · Tableaux en mémoire¶
Un tableau NumPy est un bloc de nombres contigus¶
Contrairement à une liste Python (qui stocke des références dispersés en mémoire), un tableau NumPy est un bloc contigu (contiguous block) d'octets :
liste Python [ref→12500] [ref→13600] [ref→11900] ← pointeurs dispersés
↓ ↓ ↓
numpy array | 12500 | 13600 | 11900 | ← octets côte à côte
adresse base
C'est pourquoi NumPy est rapide : le processeur peut charger plusieurs valeurs en une seule opération mémoire (vectorisation matérielle).
Conséquence pratique : un tableau NumPy a une taille fixe à la création.
Ajouter un élément (np.append) crée un nouveau tableau — c'est coûteux. Préféreznp.zeros(n)puis remplissez.
dtype — précision et empreinte mémoire¶
On peut être plus précis que float, int etc pour nos dtype. On peut spécifier aussi le niveau de précision :
dtype |
Taille | Plage / Précision | Usage |
|---|---|---|---|
float64 |
8 octets | ~15 chiffres significatifs | défaut — calcul scientifique |
float32 |
4 octets | ~7 chiffres significatifs | GPU, grandes simulations |
int64 |
8 octets | −9.2×10¹⁸ à +9.2×10¹⁸ | indices, compteurs |
int32 |
4 octets | −2.1×10⁹ à +2.1×10⁹ | |
bool |
1 octet | True / False |
masques |
a = np.array([1.0, 2.0, 3.0]) # float64 par défaut
b = np.array([1.0, 2.0, 3.0], dtype=np.float32)
c = np.zeros(5, dtype=np.int32)
print(a.dtype) # float64
print(a.nbytes) # 24 octets (3 × 8)
print(b.nbytes) # 12 octets (3 × 4)
⚠️ Convertir
float64→float32peut introduire des erreurs d'arrondi silencieuses.
Pour le travail numérique ordinaire, restez enfloat64.
Ordre de stockage : ligne-majeur (row-major) C vs colonne-majeur (column-major) Fortran¶
Un tableau 2D est stocké en mémoire comme une séquence 1D. L'ordre dépend de la convention :
A = np.array([[1, 2, 3],
[4, 5, 6]])
Ordre C (ligne-majeur) — par défaut dans NumPy :
1 2 3 4 5 6 — les colonnes de la même ligne sont adjacentes.
Ordre F (colonne-majeur) — Fortran, MATLAB :
1 4 2 5 3 6 — les lignes de la même colonne sont adjacentes.
A_c = np.array([[1,2,3],[4,5,6]], order='C') # order='C' ligne-majeur
A_f = np.array([[1,2,3],[4,5,6]], order='F') # order='F' colonne-majeur
A_c.flags['C_CONTIGUOUS'] # → True
A_f.flags['F_CONTIGUOUS'] # → True
Impact sur la performance :
En ordre C, parcourir ligne par ligne est rapide (cache-friendly) ; parcourir colonne par colonne est lent.
La règle pratique : effectuez les opérations dans le sens des lignes (axis=1) sur des tableaux C-order.
Foulées (strides) — comment NumPy navigue en mémoire¶
Les foulées (strides) décrivent le nombre d'octets à avancer pour passer d'un élément au suivant dans chaque dimension :
A = np.array([[1.0, 2.0, 3.0],
[4.0, 5.0, 6.0]]) # dtype=float64 → 8 octets par élément
A.strides # → (24, 8)
# 24 = 3 colonnes × 8 octets = avancer d'une ligne
# 8 = 1 élément × 8 octets = avancer d'une colonne
Découpage et transposition ne copient pas la mémoire — ils changent simplement les foulées :
ligne0 = A[0, :] # foulées (8,) — vue sur les 3 premiers éléments
col1 = A[:, 1] # foulées (24,) — saute de 24 octets entre éléments
AT = A.T # foulées deviennent (8, 24) — transposée sans copie
np.shares_memory(A, AT) # → True — même mémoire !
np.shares_memory(a, b)— (partage de mémoire) — retourneTruesiaetbpointent vers la même mémoire.
Vues vs copies — np.shares_memory()¶
Vue (view) : un tableau qui pointe vers la même mémoire qu'un autre.
Modifier la vue modifie l'original — et vice versa.
Copie (copy) : un nouveau tableau indépendant.
| Opération | Résultat |
|---|---|
a[2:5] |
vue |
a[::2] |
vue |
a.T |
vue |
a.reshape(n, m) |
vue si possible, copie sinon |
a[masque_booleen] |
copie — indexation "fancy" (faintaisie) |
a[[0, 2, 4]] |
copie — indexation par liste |
a.copy() |
copie explicite |
a = np.array([10, 20, 30, 40, 50])
v = a[1:4] # vue
v[0] = 999
print(a) # → [ 10 999 30 40 50] — original modifié !
c = a[1:4].copy() # copie explicite
c[0] = 0
print(a) # → [ 10 999 30 40 50] — original intact
Règle : si vous voulez modifier un sous-tableau sans affecter l'original, utilisez
.copy().
# ── Démonstration : vues, copies & strides ────────────────────────
import numpy as np
A = np.array([[1.0, 2.0, 3.0],
[4.0, 5.0, 6.0],
[7.0, 8.0, 9.0]])
print("A.strides :", A.strides) # (24, 8)
# Vue via découpage
col1 = A[:, 1]
print("col1 :", col1)
print("Partage mémoire A/col1 :", np.shares_memory(A, col1))
col1[0] = 99
print("A après modification via col1 :", A)
# Transposée — vue
AT = A.T
print("Partage mémoire A/AT :", np.shares_memory(A, AT))
# Copie explicite — indépendante
col1_copie = A[:, 1].copy()
col1_copie[0] = 0
print("A après modification via copie :", A) # inchangé
# Indexation par masque booléen — toujours une copie
sel = A[A > 5]
print("Partage mémoire A/sel (masque) :", np.shares_memory(A, sel))
⚠️ Pièges liés aux vues¶
1. Modifier une colonne extraite sans .copy()
data = np.array([[500, 22, 245],
[510, 20, 251],
[495, 24, 238]])
durete = data[:, 2] # vue sur la colonne 2
durete[0] = 0 # MODIFIE data[0, 2] !
print(data[0]) # → [500 22 0] — surprise !
durete = data[:, 2].copy() # copie indépendante — sans risque
2. reshape peut retourner une vue ou une copie selon les foulées
a = np.arange(12)
b = a.reshape(3, 4) # vue (mémoire contiguë)
c = a[::2].reshape(3, 2) # copie (foulées non contiguës après ::2)
np.shares_memory(a, b) # True
np.shares_memory(a, c) # False
3. .T sur un tableau 1D ne fait rien
v = np.array([1, 2, 3]) # shape (3,)
v.T # shape (3,) — identique ! Un vecteur n'a pas de "ligne" vs "colonne"
v.reshape(1, 3) # shape (1, 3) — vecteur ligne
v.reshape(3, 1) # shape (3, 1) — vecteur colonne
À vous de faire les exercices 1.1 et 1.2 !¶
2 · Construction et opérations matricielles¶
Construire des matrices utiles — np.eye(), np.diag()¶
import numpy as np
np.eye(3) # eye = identité → matrice identité 3×3 (float64)
# [[1. 0. 0.]
# [0. 1. 0.]
# [0. 0. 1.]]
np.diag([4, 5, 6]) # diag = diagonal → matrice diagonale à partir d'un vecteur
# [[4 0 0]
# [0 5 0]
# [0 0 6]]
np.diag(A) # Si A est 2D, extrait la diagonale
# Changer la forme — reshape
a = np.arange(6) # [0 1 2 3 4 5] shape (6,)
A = a.reshape(2, 3) # reshape = reformer [[0 1 2][3 4 5]] shape (2, 3)
A = a.reshape(3, -1) # -1 = NumPy calcule la taille automatiquement → (3, 2)
Assembler des tableaux — np.stack(), np.block()¶
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
np.stack([a, b], axis=0) # stack = empiler → shape (2, 3)
# [[1 2 3]
# [4 5 6]]
np.stack([a, b], axis=1) # → shape (3, 2)
# [[1 4]
# [2 5]
# [3 6]]
np.block() (bloc) assemble des sous-matrices comme des briques :
A11 = np.eye(2)
A12 = np.zeros((2, 2))
A21 = np.zeros((2, 2))
A22 = np.ones((2, 2))
M = np.block([[A11, A12],
[A21, A22]])
# Matrice 4×4 par blocs
np.blockest particulièrement utile pour construire des matrices de rigidité assemblées en éléments finis.
Transposée, trace & norme¶
A = np.array([[1, 2, 3],
[4, 5, 6]]) # shape (2, 3)
A.T # transposée — shape (3, 2), c'est une vue
np.trace(A) # trace → somme de la diagonale principale (nécessite tableau carré)
np.linalg.norm(A) # norme de Frobenius : sqrt(sum(aij²)) — norme par défaut
np.linalg.norm(A, ord=1) # norme-1 : max des sommes de colonnes
np.linalg.norm(A, ord=np.inf) # norme-inf : max des sommes de lignes
Pour les vecteurs :
v = np.array([3.0, 4.0])
np.linalg.norm(v) # → 5.0 (norme euclidienne — Pythagore)
np.linalg.norm(v, ord=1) # → 7.0 (norme-1 : somme des valeurs absolues)
⚠️ Multiplication matricielle : @ vs *¶
C'est la distinction la plus importante de cette séance :
| Opérateur | Signification | Exemple |
|---|---|---|
A * B |
multiplication élément par élément | [1,2] * [3,4] → [3, 8] |
A @ B |
multiplication matricielle (produit de matrices) | [1, 2] @ [3, 4].T → 11 |
A = np.array([[1, 2],
[3, 4]])
B = np.array([[5, 6],
[7, 8]])
A * B # élément par élément : [[1*5, 2*6], [3*7, 4*8]] → [[5,12],[21,32]]
A @ B # produit matriciel : [[1*5+2*7, 1*6+2*8], [3*5+4*7, 3*6+4*8]] → [[19,22],[43,50]]
Pour un vecteur u (shape (n,)) et une matrice K (shape (n,n)) :
F = K @ u # produit matrice-vecteur → shape (n,)
Erreur classique : utiliser
*là où@est attendu — Python ne lève aucune erreur si les formes sont identiques, mais le résultat est silencieusement faux.
# ── Démonstration : opérations matricielles ────────────────────
import numpy as np
# Matrice de rigidité d'un système à 3 ressorts en série
# k1=100, k2=150, k3=80 N/mm
k1, k2, k3 = 100, 150, 80
# Matrice de rigidité assemblée (tridiagonale)
K = np.array([[ k1, -k1, 0 ],
[-k1, k1+k2, -k2 ],
[ 0, -k2, k2+k3]])
print("K =", K)
# Vecteur déplacement (mm)
u = np.array([0.5, 0.3, 0.1])
# Vecteur force = K @ u (N)
F = K @ u
print("Forces F = K @ u :", F, "N")
# Norme du vecteur force
print("||F|| :", round(np.linalg.norm(F), 4), "N")
# Comparer @ et *
print("K @ u :", K @ u)
print("K * u :", K * u, " ← FAUX pour la mécanique !")
À vous de faire les exercices 2.1 et 2.2 !¶
☕ Pause — 10 minutes¶
3 · Précision numérique & virgule flottante¶
Pourquoi 0.1 + 0.2 ≠ 0.3 ?¶
0.1 + 0.2 == 0.3 # → False ← pas un bug Python, c'est de la physique informatique
0.1 + 0.2 # → 0.30000000000000004
Raison : la plupart des fractions décimales n'ont pas de représentation exacte en binaire — exactement comme 1/3 n'a pas de représentation exacte en décimal.
Les ordinateurs travaillent en base 2, et des nombres comme 0.1 = 1/10 sont des fractions périodiques en base 2 :
$$0{,}1_{10} = 0{,}0001100110011\ldots_2$$
Ce comportement est universel — il existe dans tous les langages qui utilisent la norme IEEE 754 (i.e. en pratique, tous les langages).
IEEE 754 — représentation interne (brief)¶
Un float64 (64 bits) est encodé en trois champs :
[signe : 1 bit][exposant : 11 bits][mantisse : 52 bits]
- Signe : 0 = positif, 1 = négatif
- Exposant : puissance de 2 (bias 1023)
- Mantisse : partie fractionnaire en base 2
Valeur représentée : $(-1)^s \times 1{,}\text{mantisse} \times 2^{e-1023}$
Ce qu'il faut retenir :
float64a environ 15–16 chiffres significatifs en base 10- Il y a un nombre fini de valeurs représentables — entre deux valeurs consécutives, il y a toujours un écart
- Cet écart minimal s'appelle l'epsilon machine (machine epsilon)
import numpy as np
eps = np.finfo(np.float64).eps # finfo = float info
print(eps) # → 2.220446049250313e-16
# Signification : le plus petit δ tel que 1.0 + δ ≠ 1.0 en float64
Annulation catastrophique (catastrophic cancellation)¶
L'annulation catastrophique survient lorsqu'on soustrait deux nombres presque égaux — les chiffres significatifs s'annulent et le résultat perd de la précision.
Exemple :
a = 1.000000001
b = 1.000000000
c = a - b # mathématiquement : 1e-9
# résultat : 1.0000000009313226e-09 ← erreur relative ~7 %
Si vous calculez avec des points proches, l'erreur peut être hors tolérances.
Règle pratique : quand les deux termes d'une soustraction sont proches, cherchez une formulation algébriquement équivalente qui évite la soustraction.
# Mauvais : (√(x+h) - √x) / h pour h petit → annulation catastrophique
# Meilleur : 1 / (√(x+h) + √x) — formule conjuguée, algébriquement identique
# ── Démonstration : précision virgule flottante ───────────────
import numpy as np
# 0.1 + 0.2 ≠ 0.3
print("0.1 + 0.2 == 0.3 :", 0.1 + 0.2 == 0.3)
print("0.1 + 0.2 :", repr(0.1 + 0.2))
# Epsilon machine
eps = np.finfo(np.float64).eps
print("Epsilon machine float64 :", eps)
print("1.0 + eps/2 == 1.0 :", 1.0 + eps/2 == 1.0) # True — trop petit pour être vu
print("1.0 + eps == 1.0 :", 1.0 + eps == 1.0) # False
# (1 + petit) - 1
print("(1 + 1e-16) - 1 :", (1 + 1e-16) - 1) # → 0.0 perte totale !
print("(1 + 1e-8 ) - 1 :", (1 + 1e-8) - 1) # → 9.999999994515327e-09
# Somme 1/n² — avant vs arrière
N = 100000
n = np.arange(1, N + 1)
avant = np.sum(1.0 / n**2) # grand → petit : moins précis
arriere = np.sum(1.0 / n[::-1]**2) # petit → grand : plus précis
vrai = np.pi**2 / 6
print("Somme 1/n² (avant) :", avant, " erreur :", abs(avant - vrai))
print("Somme 1/n² (arrière) :", arriere, " erreur :", abs(arriere - vrai))
print("π²/6 (scipy) :", vrai)
⚠️ Ne jamais comparer des flottants avec ==¶
# MAUVAIS — toujours risqué avec des flottants
if sigma == 500.0:
print("exactement 500 MPa") # peut ne jamais s'exécuter à cause des erreurs d'arrondi
# BON — utiliser np.isclose() ou np.allclose()
np.isclose(a, b, rtol=1e-5, atol=1e-8)
# isclose = proche de → True si |a - b| ≤ atol + rtol * |b|
np.allclose(arr1, arr2) # allclose = tous proches → True si tous les éléments sont proches
| Fonction | Usage |
|---|---|
np.isclose(a, b) |
comparer deux scalaires ou deux tableaux élément par élément |
np.allclose(A, B) |
vérifier qu'un tableau entier est "égal" à un autre |
# Exemple typique : vérifier qu'une solution est correcte
x = np.linalg.solve(A, b)
np.allclose(A @ x, b) # True si la solution est bonne — la bonne façon de vérifier
À vous de faire les exercices 3.1 et 3.2 !¶
4 · Analyse d'erreur & propagation d'incertitude¶
Erreur absolue et erreur relative¶
$$\text{erreur absolue} = |x_{\text{calculé}} - x_{\text{vrai}}|$$
$$\text{erreur relative} = \frac{|x_{\text{calculé}} - x_{\text{vrai}}|}{|x_{\text{vrai}}|}$$
Quand utiliser laquelle ?
| Situation | Métrique appropriée |
|---|---|
| Comparer des valeurs de même grandeur | Erreur absolue |
| Comparer des valeurs d'ordre de grandeur très différents | Erreur relative |
| Spécifier la qualité d'une mesure | Erreur relative (%) — ex. "précision de 0.1%" |
x_vrai = 210000.0 # GPa — module de Young acier
x_calc = 209876.5
err_abs = abs(x_calc - x_vrai) # → 123.5
err_rel = abs(x_calc - x_vrai) / x_vrai # → 5.88e-4 soit 0.059 %
Propagation d'incertitude¶
Si une grandeur calculée $y = f(x_1, x_2, \ldots)$ dépend de mesures avec des incertitudes $\delta x_i$, l'incertitude sur $y$ est :
$$\delta y = \sqrt{\left(\frac{\partial f}{\partial x_1}\delta x_1\right)^2 + \left(\frac{\partial f}{\partial x_2}\delta x_2\right)^2 + \cdots}$$
Exemple — contrainte d'ingénierie : $\sigma = F / A = F / (\pi (d/2)^2)$
$$\frac{\delta\sigma}{\sigma} = \sqrt{\left(\frac{\delta F}{F}\right)^2 + \left(\frac{2\,\delta d}{d}\right)^2}$$
F, delta_F = 15000, 50 # N — force ± 50 N
d, delta_d = 10.0, 0.01 # mm — diamètre ± 0.01 mm
import numpy as np
A = np.pi * (d / 2)**2
sigma = F / A
delta_sigma = sigma * np.sqrt((delta_F/F)**2 + (2*delta_d/d)**2)
print("σ =", round(sigma, 2), "±", round(delta_sigma, 2), "MPa")
Règle des chiffres significatifs : on ne reporte jamais plus de chiffres significatifs que ce que justifie l'incertitude.
Si $\delta\sigma \approx 5$ MPa, reporter $\sigma = 190$ MPa (3 chiffres), pas $\sigma = 190{,}99$ MPa.
# ── Démonstration : propagation d'incertitude ─────────────────
import numpy as np
def contrainte_avec_incertitude(F, d, delta_F, delta_d):
"""
Calcule la contrainte d'ingénierie et son incertitude propagée.
Paramètres
----------
F : float Force (N)
d : float Diamètre (mm)
delta_F : float Incertitude sur F (N)
delta_d : float Incertitude sur d (mm)
Retourne
--------
sigma : float Contrainte (MPa)
delta_sigma : float Incertitude sur la contrainte (MPa)
"""
A = np.pi * (d / 2)**2
sigma = F / A
delta_sigma = sigma * np.sqrt((delta_F / F)**2 + (2 * delta_d / d)**2)
return sigma, delta_sigma
# Données
forces = np.array([15000, 16200, 14800, 15500]) # N
diametr = np.array([10.00, 9.98, 10.02, 10.01]) # mm
dF = 50 # N — incertitude sur la force
dd = 0.01 # mm — incertitude sur le diamètre
for i in range(len(forces)):
s, ds = contrainte_avec_incertitude(forces[i], diametr[i], dF, dd)
print("Eprouvette", i + 1, ": sigma =", round(s, 1), "+/-", round(ds, 1), "MPa")
# Vérification vectorisée
A_vec = np.pi * (diametr / 2)**2
sig_vec = forces / A_vec
ds_vec = sig_vec * np.sqrt((dF / forces)**2 + (2 * dd / diametr)**2)
print("Vectorisé :")
print("sigma :", np.round(sig_vec, 1), "MPa")
print("delta :", np.round(ds_vec, 2), "MPa")
À vous de faire les exercices 4.1 et 4.2 !¶
Points clés de la séance¶
| Concept | Ce qu'il faut retenir |
|---|---|
| Mémoire contiguë | Un tableau NumPy est un bloc d'octets — c'est ce qui le rend rapide |
dtype |
Détermine précision et mémoire ; rester en float64 par défaut |
| Ordre C vs F | NumPy est C-order : itérer sur les lignes est plus efficace |
| Foulées | .T, a[::2], a[:, 1] sont des vues — pas de copie |
| Vues vs copies | Modifier une vue modifie l'original ; .copy() pour l'indépendance |
@ vs * |
@ = produit matriciel ; * = élément par élément — ne pas confondre |
0.1 + 0.2 ≠ 0.3 |
Pas un bug Python — c'est la norme IEEE 754 |
| Epsilon machine | np.finfo(np.float64).eps ≈ 2.2×10⁻¹⁶ |
| Annulation catastrophique | Soustraction de valeurs proches → perte de précision |
np.isclose |
Toujours comparer des flottants avec isclose, jamais == |
| Propagation d'incertitude | $\delta y = \sqrt{\sum(\partial f / \partial x_i \cdot \delta x_i)^2}$ |