import os
import numpy as np
import matplotlib.pyplot as plt

# SOLUTION
def estimer_pi(n, seed=0):
    rng = np.random.default_rng(seed)
    points = rng.uniform(-1, 1, size=(n, 2))
    dedans = (points**2).sum(axis=1) <= 1
    return 4 * dedans.mean()

estimer_pi(1_000_000)

np.float64(3.141576)

assert 3.0 < estimer_pi(1000) < 3.3
assert abs(estimer_pi(1_000_000, seed=0) - np.pi) < 0.01
assert estimer_pi(5000, seed=42) == estimer_pi(5000, seed=42)
assert estimer_pi(5000, seed=1) != estimer_pi(5000, seed=2)
print("ok")

ok

# SOLUTION
rng = np.random.default_rng(0)
pts = rng.uniform(-1, 1, size=(3000, 2))
masque = (pts**2).sum(axis=1) <= 1
plt.figure(figsize=(5, 5))
plt.scatter(pts[masque, 0], pts[masque, 1], s=4, label="dans le disque")
plt.scatter(pts[~masque, 0], pts[~masque, 1], s=4, label="dehors")
plt.axis("equal")
plt.legend(loc="upper right")
plt.title(f"3000 points : π ≈ {4 * masque.mean():.4f}")
plt.show()

# SOLUTION
ns = np.logspace(2, 6, 9).astype(int)
graines = range(5)
erreurs = np.array([[abs(estimer_pi(n, seed=s) - np.pi) for s in graines] for n in ns]).mean(axis=1)

plt.loglog(ns, erreurs, "o-", label="erreur moyenne sur 5 graines")
plt.loglog(ns, 1 / np.sqrt(ns), "k--", label=r"$1/\sqrt{n}$")
plt.xlabel("nombre de points n")
plt.ylabel(r"$|\hat\pi_n - \pi|$")
plt.title("Convergence de l'estimation Monte-Carlo de π")
plt.legend()
plt.grid(True, which="both", alpha=0.3)
plt.savefig("pi.png", dpi=120)
plt.show()

assert ns.shape == (9,) and ns[0] == 100 and ns[-1] == 1_000_000
assert erreurs.shape == (9,)
assert erreurs[-1] < erreurs[0]
assert erreurs[-1] < 0.01
assert os.path.exists("pi.png")
print("ok")

ok

# SOLUTION
K, n = 500, 1000
rng = np.random.default_rng(0)
tirages = rng.integers(0, 4, size=(K, n))
directions = np.array([[1, 0], [-1, 0], [0, 1], [0, -1]])
pas = directions[tirages]
positions = pas.cumsum(axis=1)
print(pas.shape, positions.shape)

(500, 1000, 2) (500, 1000, 2)

assert tirages.shape == (500, 1000)
assert directions.shape == (4, 2)
assert pas.shape == (500, 1000, 2)
assert positions.shape == (500, 1000, 2)
assert np.all(np.abs(pas).sum(axis=2) == 1)          # chaque pas a une longueur 1
assert np.array_equal(positions[:, 0, :], pas[:, 0, :])
assert np.array_equal(positions[:, -1, :], pas.sum(axis=1))
print("ok")

ok

# SOLUTION
finales = positions[:, -1, :]
distances = np.sqrt((finales**2).sum(axis=1))
print("RMS :", np.sqrt((distances**2).mean()), "  sqrt(n) =", np.sqrt(n))
print("moyenne :", distances.mean(), "  0.886 sqrt(n) =", np.sqrt(np.pi * n) / 2)
print("proportion à moins de 20 :", (distances < 20).mean())
dist_t = np.sqrt((positions**2).sum(axis=2))
print("distance max atteinte (tous marcheurs, tous instants) :", dist_t.max())
print("distance max atteinte, marcheur par marcheur (5 premiers) :", dist_t.max(axis=1)[:5])

RMS : 31.73824191728332   sqrt(n) = 31.622776601683793
moyenne : 28.11196907157812   0.886 sqrt(n) = 28.024956081989643
proportion à moins de 20 : 0.324
distance max atteinte (tous marcheurs, tous instants) : 86.03487664894976
distance max atteinte, marcheur par marcheur (5 premiers) : [21.21320344 45.79301257 44.68780594 28.0713377  19.84943324]

assert finales.shape == (500, 2) and distances.shape == (500,)
rms = np.sqrt((distances**2).mean())
assert 0.85 * np.sqrt(n) < rms < 1.15 * np.sqrt(n)
assert dist_t.shape == (500, 1000)
assert np.all(dist_t.max(axis=1) >= distances)
print("ok")

ok

# SOLUTION
plt.figure(figsize=(6, 6))
for k in range(5):
    plt.plot(positions[k, :, 0], positions[k, :, 1], lw=0.8, label=f"marcheur {k}")
plt.plot(0, 0, "ko", ms=8, label="origine")
plt.axis("equal")
plt.legend()
plt.title(f"5 marches aléatoires de {n} pas")
plt.savefig("marche.png", dpi=120)
plt.show()

t = np.arange(1, n + 1)
plt.plot(t, dist_t.mean(axis=0), label="distance moyenne (500 marcheurs)")
plt.plot(t, 0.886 * np.sqrt(t), "k--", label=r"$0{,}886\sqrt{t}$")
plt.xlabel("pas t"); plt.ylabel("distance à l'origine")
plt.legend(); plt.title("Éloignement moyen au cours du temps")
plt.show()

assert os.path.exists("marche.png")
print("ok")

ok

# SOLUTION
retour = (positions == 0).all(axis=2).any(axis=1)
print("proportion de retours à l'origine :", retour.mean())

proportion de retours à l'origine : 0.7

assert retour.shape == (500,) and retour.dtype == bool
assert 0.5 < retour.mean() < 0.95
print("ok")

ok

# SOLUTION
def grille(xmin, xmax, ymin, ymax, nx, ny):
    x = np.linspace(xmin, xmax, nx)
    y = np.linspace(ymin, ymax, ny)
    return x[np.newaxis, :] + 1j * y[:, np.newaxis]

grille(-1, 1, -1, 1, 3, 2)

array([[-1.-1.j,  0.-1.j,  1.-1.j],
       [-1.+1.j,  0.+1.j,  1.+1.j]])

g = grille(-2, 1, -1.5, 1.5, 4, 3)
assert g.shape == (3, 4)
assert g.dtype == complex
assert g[0, 0] == -2 - 1.5j and g[2, 3] == 1 + 1.5j and g[1, 0] == -2 + 0j
print("ok")

ok

# SOLUTION
def mandelbrot(c, max_iter=100):
    z = np.zeros_like(c)
    n_iter = np.zeros(c.shape, dtype=int)
    for k in range(max_iter):
        z = z**2 + c
        encore = np.abs(z) <= 2
        n_iter += encore
    return n_iter

with np.errstate(over="ignore", invalid="ignore"):
    m = mandelbrot(grille(-2, 1, -1.5, 1.5, 300, 300), max_iter=50)
m.shape, m.min(), m.max()

((300, 300), np.int64(0), np.int64(50))

with np.errstate(over="ignore", invalid="ignore"):
    m = mandelbrot(grille(-2, 1, -1.5, 1.5, 300, 300), max_iter=50)
assert m.shape == (300, 300)
assert m.dtype.kind in "iu"
assert m.max() == 50 and m.min() >= 0
assert m[150, 100] == 50          # c ≈ -1 : dans l'ensemble
assert m[0, 0] <= 2               # c = -2 - 1.5i : diverge tout de suite
assert m[150, 250] <= 5           # c ≈ 0.5 : diverge vite
assert 0.14 < (m == 50).mean() < 0.21   # aire de l'ensemble ≈ 1.5 sur un cadre de 9
print("ok")

ok

# SOLUTION
xmin, xmax, ymin, ymax = -2, 1, -1.5, 1.5
with np.errstate(over="ignore", invalid="ignore"):
    M = mandelbrot(grille(xmin, xmax, ymin, ymax, 900, 900), max_iter=100)

plt.figure(figsize=(7, 7))
plt.imshow(M, extent=[xmin, xmax, ymin, ymax], origin="lower", cmap="magma")
plt.colorbar(label="itérations", shrink=0.8)
plt.xlabel("Re(c)"); plt.ylabel("Im(c)")
plt.title("Ensemble de Mandelbrot (max_iter = 100)")
plt.savefig("mandelbrot.png", dpi=120)
plt.show()

assert M.shape == (900, 900)
assert os.path.exists("mandelbrot.png")
print("ok")

ok

# SOLUTION
def mandelbrot_masque(c, max_iter=100):
    z = np.zeros_like(c)
    n_iter = np.zeros(c.shape, dtype=int)
    actif = np.ones(c.shape, dtype=bool)
    for k in range(max_iter):
        z[actif] = z[actif]**2 + c[actif]
        encore = np.abs(z[actif]) <= 2   # parmi les points actifs, ceux qui restent bornés
        n_iter[actif] += encore
        actif[actif] = encore            # on désactive les autres
    return n_iter

c2 = grille(-2, 1, -1.5, 1.5, 300, 300)
with np.errstate(over="ignore", invalid="ignore"):
    ref = mandelbrot(c2, 50)
print(np.array_equal(mandelbrot_masque(c2, 50), ref))

zx = (-0.8, -0.7, 0.05, 0.15)
Z = mandelbrot_masque(grille(*zx, 600, 600), max_iter=300)
plt.figure(figsize=(7, 7))
plt.imshow(Z, extent=zx, origin="lower", cmap="twilight")
plt.colorbar(label="itérations", shrink=0.8)
plt.xlabel("Re(c)"); plt.ylabel("Im(c)")
plt.title("Zoom : vallée des hippocampes (max_iter = 300)")
plt.show()

True

# CORRIGE-ONLY
with np.errstate(over="ignore", invalid="ignore"):
    %timeit mandelbrot(c2, 50)
%timeit mandelbrot_masque(c2, 50)

8.65 ms ± 163 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

10.5 ms ± 74.8 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Partie	Contenu	Points
A	Monte-Carlo : estimation de π, convergence, `pi.png`	6
B	Marches aléatoires 2D vectorisées, `marche.png`	6
C	Ensemble de Mandelbrot, `mandelbrot.png`	6
–	Notebook propre : exécuté sans erreur, code vectorisé, figures légendées, réponses rédigées	2
Total		20
Bonus	Zoom, palette, itération sans avertissement	+2

D2 – Monte-Carlo & Mandelbrot¶

Barème¶

Partie A – Estimer π par la méthode de Monte-Carlo (6 pts)¶

Partie B – Marches aléatoires dans le plan (6 pts)¶

Partie C – L'ensemble de Mandelbrot (6 pts)¶

Bonus (+2)¶