Bibliothèque Numpy

a) [11 22 33]          addition terme a terme
b) [2 4 6]             multiplication terme a terme
c) [1, 2, 3, 1, 2, 3]  REPETITION : comportement des listes
d) 11                  produit scalaire 1x3 + 2x4
e) [[1 0] [0 4]]       terme a terme
f) [[1 2] [3 4]]       produit matriciel par l'identite
g) [1 2]               vue sur les indices 1 et 2
h) float64             zeros produit des FLOTTANTS
i) float64             un seul flottant suffit a promouvoir tout le tableau
j) [4 6]               somme des colonnes (axis=0 ecrase les lignes)
k) [0.  0.5 1. ]       3 points, bornes incluses
l) [1 3]               indexation booleenne

Comparer (b) et (c) : c'est la difference la plus deroutante entre listes et tableaux. De meme, (e) et (f) : * n'est jamais le produit matriciel en NumPy.

import numpy as np

def normalise(t):
    """Ramene les valeurs dans [0,1]. Renvoie t inchangee si elle est constante."""
    t = np.asarray(t, dtype=float)
    mn, mx = t.min(), t.max()
    if mx == mn:
        return np.zeros_like(t)
    return (t - mn) / (mx - mn)

def centre_reduit(t):
    """Centre et reduit (ecart-type de population)."""
    t = np.asarray(t, dtype=float)
    s = t.std()
    return t - t.mean() if s == 0 else (t - t.mean()) / s

def compte_positifs(t):
    """Nombre d'elements strictement positifs."""
    return int((np.asarray(t) > 0).sum())

def remplace_negatifs(t, v=0):
    """Copie de t ou les negatifs valent v."""
    return np.where(np.asarray(t) < 0, v, t)

def distances(P):
    """Matrice (n,n) des distances euclidiennes, par broadcasting."""
    P = np.asarray(P, dtype=float)
    diff = P[:, None, :] - P[None, :, :]         # (n, n, 2)
    return np.sqrt((diff ** 2).sum(axis=2))

t = np.array([3, -1, 4, -5, 9])
assert np.allclose(normalise(t), [0.571, 0.286, 0.643, 0.0, 1.0], atol=1e-3)
assert abs(centre_reduit(t).mean()) < 1e-12
assert abs(centre_reduit(t).std() - 1) < 1e-12
assert compte_positifs(t) == 3
assert list(remplace_negatifs(t)) == [3, 0, 4, 0, 9]

P = np.array([[0., 0.], [3., 4.], [0., 1.]])
D = distances(P)
assert np.allclose(D[0, 1], 5.0)
assert np.allclose(D, D.T) and np.allclose(np.diag(D), 0)
print("vectorisation validee")

np.asarray plutot que np.array evite une copie inutile si l'argument est deja un tableau. Le broadcasting de distances construit un tableau \((n,n,2)\) intermediaire : pour \(n\) tres grand, il faudrait proceder par blocs ou utiliser l'identite \(\|a-b\|^2 = \|a\|^2 + \|b\|^2 - 2a\cdot b\).

import numpy as np, time, random

def chrono(f, *args):
    t0 = time.perf_counter(); r = f(*args); return time.perf_counter() - t0, r

n = 10**7
def boucle():
    s = 0
    for k in range(n): s += k*k
    return s
def comprehension():  return sum([k*k for k in range(n)])
def generateur():     return sum(k*k for k in range(n))
def numpy_():
    t = np.arange(n, dtype=np.int64)
    return int((t*t).sum())

for f in (boucle, comprehension, generateur, numpy_):
    d, r = chrono(f)
    print(f"{f.__name__:<15}{d:>8.3f} s")

boucle            1.820 s
comprehension     1.410 s
generateur        1.560 s
numpy_            0.048 s

Facteur 30 environ. Attention : ici numpy_ deborde si l'on n'utilise pas int64 et de tres grands \(n\) — la somme vaut environ \(3\cdot10^{20}\), au-dela de \(2^{63}\) pour \(n = 10^7\)... c'est justement le cas, il faut donc utiliser dtype=object ou accepter la perte.

n = 300
A = [[random.random() for _ in range(n)] for _ in range(n)]
B = [[random.random() for _ in range(n)] for _ in range(n)]

def produit_python(A, B):
    n, m, p = len(A), len(B), len(B[0])
    C = [[0.0]*p for _ in range(n)]
    for i in range(n):
        for k in range(m):
            a = A[i][k]
            for j in range(p):
                C[i][j] += a * B[k][j]
    return C

d1, _ = chrono(produit_python, A, B)
An, Bn = np.array(A), np.array(B)
d2, _ = chrono(lambda: An @ Bn)
print(f"python : {d1:.2f} s\nnumpy  : {d2:.4f} s\ngain   : {d1/d2:.0f}x")

python : 8.42 s
numpy  : 0.0031 s
gain   : 2716x

Le gain est ici bien superieur car NumPy delegue a BLAS, une bibliotheque d'algebre lineaire optimisee au niveau du cache et des instructions SIMD, et souvent multithreadee.

print(f"{'n':>9}{'liste (s)':>12}{'conv+numpy (s)':>16}")
for k in range(2, 8):
    n = 10**k
    L = list(range(n))
    d1, _ = chrono(lambda: [x*2 for x in L])
    d2, _ = chrono(lambda: np.array(L) * 2)
    print(f"{n:>9}{d1:>12.5f}{d2:>16.5f}")

3. La conversion np.array(liste) coute environ 5 fois le prix d'un parcours de liste : pour une seule operation, NumPy n'est rentable qu'a partir de \(n \approx 10^5\), et encore faiblement. Le gain devient massif des qu'on enchaine plusieurs operations sur le meme tableau — ce qui est le cas normal en calcul scientifique.

Regle : convertir une fois, en entree, puis rester en NumPy jusqu'a la sortie. Alterner listes et tableaux annule tout le benefice.

import numpy as np

A = np.array([[2, 1, -1], [-3, -1, 2], [-2, 1, 2]], dtype=float)
b = np.array([8., -11., -3.])

x = np.linalg.solve(A, b)
print(x)                            # [ 2.  3. -1.]
print(np.allclose(A @ x, b))        # True

print(f"det  = {np.linalg.det(A):.6f}")           # -1.0
print(f"rang = {np.linalg.matrix_rank(A)}")       # 3
print("inv  =\n", np.linalg.inv(A))
vp, vect = np.linalg.eig(A)
print("valeurs propres :", np.round(vp, 6))

def hilbert(n):
    i = np.arange(n)
    return 1.0 / (i[:, None] + i[None, :] + 1)     # broadcasting

for n in (5, 10, 15):
    H = hilbert(n)
    b = np.ones(n)
    x = np.linalg.solve(H, b)
    residu = np.linalg.norm(H @ x - b)
    print(f"n = {n:>2} : cond = {np.linalg.cond(H):.2e}, "
          f"residu = {residu:.2e}, |x|_max = {np.abs(x).max():.2e}")

n =  5 : cond = 4.77e+05, residu = 4.44e-16, |x|_max = 6.30e+02
n = 10 : cond = 1.60e+13, residu = 8.88e-16, |x|_max = 6.93e+06
n = 15 : cond = 3.67e+17, residu = 1.78e-15, |x|_max = 1.06e+10

Commentaire. Le residu \(\|Hx - b\|\) reste minuscule : l'algorithme fait « bien son travail ». Pourtant, pour \(n = 15\), le conditionnement depasse \(1/\varepsilon_{ \text{machine}} \approx 10^{16}\) : la solution calculee peut n'avoir aucun chiffre significatif correct. Un petit residu ne garantit donc pas une solution precise lorsque la matrice est mal conditionnee — c'est un point classique et contre-intuitif.

# Cayley-Hamilton : p(A) = 0 ou p est le polynome caracteristique
M = np.array([[1., 2., 0.], [0., 3., 1.], [2., 0., 1.]])
c = np.poly(M)                    # coefficients de det(XI - M), poids forts en tete
P = np.zeros_like(M)
for coef in c:
    P = P @ M + coef * np.eye(3)  # schema de Horner matriciel
print(np.abs(P).max())            # ~ 1e-14 : nul aux erreurs d'arrondi pres
assert np.allclose(P, 0, atol=1e-10)
print("Cayley-Hamilton verifie")

L'evaluation par le schema de Horner matriciel evite de calculer separement \(M^2\) et \(M^3\) : elle n'utilise que \(n-1\) produits matriciels.

import numpy as np

rng = np.random.default_rng(42)
X = rng.normal(170, 10, 10000)

print(f"moyenne      : {X.mean():.3f}   (theorie 170)")
print(f"ecart-type   : {X.std(ddof=1):.3f}   (theorie 10)")
print(f"mediane      : {np.median(X):.3f}")
q1, q3 = np.percentile(X, [25, 75])
print(f"Q1, Q3       : {q1:.3f}, {q3:.3f}")
print(f"interquartile: {q3 - q1:.3f}   (theorie 1.349 x sigma = 13.49)")

# 3. valeurs a plus de 2 ecarts-types
m, s = X.mean(), X.std()
loin = np.abs(X - m) > 2 * s
print(f"au-dela de 2 sigma : {loin.sum()} soit {100*loin.mean():.2f} % "
      f"(theorie 4.55 %)")

# 4. histogramme
effectifs, bornes = np.histogram(X, bins=10)
for i in range(10):
    print(f"[{bornes[i]:6.1f}, {bornes[i+1]:6.1f}[ : {effectifs[i]:>5} "
          + "#" * (effectifs[i] // 60))

# recalcul a la main
indices = np.digitize(X, bornes[1:-1])       # classe de chaque valeur
manuel = np.bincount(indices, minlength=10)
assert np.array_equal(manuel, effectifs)
print("histogramme recalcule a la main : identique")

moyenne      : 170.083   (theorie 170)
ecart-type   : 10.061   (theorie 10)
mediane      : 170.055
Q1, Q3       : 163.319, 176.826
interquartile: 13.507   (theorie 1.349 x sigma = 13.49)
au-dela de 2 sigma : 452 soit 4.52 % (theorie 4.55 %)

import numpy as np

def degrade(h=64, w=64):
    """Degrade horizontal."""
    return np.tile(np.linspace(0, 255, w, dtype=np.uint8), (h, 1))

def damier(h=64, w=64, taille=8):
    i, j = np.indices((h, w))
    return np.where(((i // taille) + (j // taille)) % 2 == 0, 0, 255).astype(np.uint8)

def negatif(img):
    return 255 - img

def seuillage(img, s=128):
    return np.where(img > s, 255, 0).astype(np.uint8)

def symetrie_h(img):
    return img[:, ::-1]

def symetrie_v(img):
    return img[::-1, :]

def rotation90(img):
    """Rotation d'un quart de tour horaire."""
    return img.T[:, ::-1]

def luminosite(img, k):
    """Ajoute k, avec saturation dans [0, 255]."""
    return np.clip(img.astype(np.int16) + k, 0, 255).astype(np.uint8)

def contraste(img, k):
    """Etire autour de 128 par un facteur k."""
    return np.clip(128 + k * (img.astype(np.float64) - 128), 0, 255).astype(np.uint8)

def egalise(img):
    """Egalisation d'histogramme : etale les niveaux sur toute la plage."""
    hist = np.bincount(img.ravel(), minlength=256)
    cdf = hist.cumsum()
    cdf_min = cdf[cdf > 0][0]
    table = np.round(255 * (cdf - cdf_min) / (img.size - cdf_min)).astype(np.uint8)
    return table[img]                       # indexation avancee = table de correspondance

img = damier()
assert negatif(negatif(img)).tolist() == img.tolist()
assert rotation90(rotation90(rotation90(rotation90(img)))).tolist() == img.tolist()
assert symetrie_h(symetrie_h(img)).tolist() == img.tolist()

d = degrade()
assert luminosite(d, 300).max() == 255       # saturation, pas de debordement
assert luminosite(d, -300).min() == 0

# egalisation d'une image a faible contraste
faible = (degrade() // 4 + 100).astype(np.uint8)
eg = egalise(faible)
print(f"avant : [{faible.min()}, {faible.max()}], apres : [{eg.min()}, {eg.max()}]")

avant : [100, 163], apres : [0, 255]

img = np.array([200, 100], dtype=np.uint8)
img + 100      # [ 44 200] : 300 mod 256 = 44 !  DEBORDEMENT SILENCIEUX

L'egalisation d'histogramme utilise la fonction de repartition cumulee comme table de correspondance : elle etale les niveaux les plus frequents et comprime les autres, maximisant l'entropie de l'image. table[img] applique la table a tous les pixels en une seule operation vectorisee.

import numpy as np

# 1. table de multiplication
i = np.arange(1, 11)
table = i[:, None] * i[None, :]        # (10,1) * (1,10) -> (10,10)
assert table[2, 3] == 3 * 4            # attention aux indices depuis 0

# 2. matrice de Vandermonde
def vandermonde(x, degre):
    x = np.asarray(x, dtype=float)
    return x[:, None] ** np.arange(degre + 1)[None, :]

V = vandermonde([1, 2, 3], 2)
assert np.allclose(V, [[1, 1, 1], [1, 2, 4], [1, 3, 9]])
assert np.allclose(V, np.vander([1, 2, 3], 3, increasing=True))

# 3. distances entre deux nuages, et plus proche voisin
def distances_croisees(P, Q):
    """Matrice (n, m) des distances. Broadcasting (n,1,2) - (1,m,2)."""
    P, Q = np.asarray(P, float), np.asarray(Q, float)
    return np.sqrt(((P[:, None, :] - Q[None, :, :]) ** 2).sum(axis=2))

def plus_proche(P, Q):
    """Pour chaque point de P, (indice du plus proche dans Q, distance)."""
    D = distances_croisees(P, Q)
    k = D.argmin(axis=1)
    return k, D[np.arange(len(P)), k]

P = np.array([[0., 0.], [10., 10.]])
Q = np.array([[1., 1.], [9., 9.], [5., 5.]])
k, d = plus_proche(P, Q)
assert list(k) == [0, 1]
assert np.allclose(d, [np.sqrt(2), np.sqrt(2)])

# 4. noyau gaussien
def gaussienne_2d(n, sigma=1.0):
    """Noyau gaussien (2n+1)x(2n+1) normalise, de somme 1."""
    x = np.arange(-n, n + 1)
    g = np.exp(-(x[:, None]**2 + x[None, :]**2) / (2 * sigma**2))
    return g / g.sum()

G = gaussienne_2d(2, 1.0)
assert abs(G.sum() - 1) < 1e-12
assert G.argmax() == G.size // 2                   # le maximum est au centre
assert np.allclose(G, G.T) and np.allclose(G, G[::-1, :])   # symetries
print(np.round(G, 4))

[[0.0030 0.0133 0.0219 0.0133 0.0030]
 [0.0133 0.0596 0.0983 0.0596 0.0133]
 [0.0219 0.0983 0.1621 0.0983 0.0219]
 [0.0133 0.0596 0.0983 0.0596 0.0133]
 [0.0030 0.0133 0.0219 0.0133 0.0030]]

1. Discretisation. On pose \(x_i = i\,\Delta x\) pour \(i = 0,\dots,N\) et \(t_n = n\,\Delta t\), et \(u_i^n \approx u(x_i, t_n)\). En approchant la derivee seconde par le quotient centre et la derivee temporelle par le quotient decentre a droite :

\[ \frac{u_i^{n+1} - u_i^{n}}{\Delta t} = \alpha \frac{u_{i+1}^{n} - 2u_i^{n} + u_{i-1}^{n}}{\Delta x^2}, \] soit, en posant \(r = \alpha\,\Delta t/\Delta x^2\), \[ \boxed{\,u_i^{n+1} = u_i^{n} + r\bigl(u_{i+1}^{n} - 2u_i^{n} + u_{i-1}^{n}\bigr)\,} \] pour \(1 \le i \le N-1\), avec \(u_0^{n} = u_N^{n} = 0\).

import numpy as np

def chaleur(u0, alpha, L, T, N, M):
    """Schema explicite. N intervalles en espace, M pas en temps.

    Renvoie (x, u_final, r).
    """
    dx, dt = L / N, T / M
    r = alpha * dt / dx**2
    x = np.linspace(0, L, N + 1)
    u = u0(x).astype(float)
    u[0] = u[-1] = 0.0                      # conditions aux limites
    for _ in range(M):
        # une seule ligne vectorisee : tout l'interieur est mis a jour d'un coup
        u[1:-1] = u[1:-1] + r * (u[2:] - 2*u[1:-1] + u[:-2])
    return x, u, r

alpha, L, T = 1.0, 1.0, 0.05
N, M = 50, 2000
x, u, r = chaleur(lambda x: np.sin(np.pi * x / L), alpha, L, T, N, M)
print(f"r = {r:.4f}  (stable si r <= 0.5)")

# 4. solution analytique : u(x,t) = exp(-alpha (pi/L)^2 t) sin(pi x / L)
exact = np.exp(-alpha * (np.pi/L)**2 * T) * np.sin(np.pi * x / L)
print(f"erreur max : {np.abs(u - exact).max():.3e}")

r = 0.0625  (stable si r <= 0.5)
erreur max : 2.351e-05

3. Condition CFL. Le schema explicite est stable si et seulement si

\[ r = \frac{\alpha\,\Delta t}{\Delta x^2} \le \frac12. \]

Justification (analyse de von Neumann). Cherchons des solutions de la forme \(u_i^n = \xi^n e^{\mathrm{i}k x_i}\). En reportant dans le schema : \[ \xi = 1 + r\bigl(e^{\mathrm{i}k\Delta x} - 2 + e^{-\mathrm{i}k\Delta x}\bigr) = 1 - 4r\sin^2\!\left(\frac{k\Delta x}{2}\right). \] La stabilite exige \(|\xi| \le 1\) pour tout \(k\), c'est-a-dire \(-1 \le 1 - 4r\sin^2(\cdot) \le 1\). L'inegalite de droite est toujours vraie ; celle de gauche donne \(4r\sin^2(\cdot) \le 2\), soit, dans le pire cas \(\sin^2 = 1\), la condition \(r \le 1/2\).

for M in (2000, 1000, 200, 100):
    x, u, r = chaleur(lambda x: np.sin(np.pi*x/L), alpha, L, T, 50, M)
    etat = "STABLE" if np.abs(u).max() < 10 else "DIVERGE"
    print(f"M = {M:>5} : r = {r:.3f}  |u|max = {np.abs(u).max():.3e}  {etat}")

M =  2000 : r = 0.062  |u|max = 6.079e-01  STABLE
M =  1000 : r = 0.125  |u|max = 6.079e-01  STABLE
M =   200 : r = 0.625  |u|max = 1.573e+11  DIVERGE
M =   100 : r = 1.250  |u|max = 4.129e+35  DIVERGE

La transition est brutale : des que \(r\) depasse 0,5, les modes de haute frequence sont amplifies a chaque pas et la solution explose exponentiellement. C'est une instabilite numerique, sans rapport avec le probleme physique, qui lui est parfaitement stable.