Bibliothèque matplotlib

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-2*np.pi, 2*np.pi, 500)

plt.figure(figsize=(9, 5))
plt.plot(x, np.sin(x),   color="#026573", linewidth=2,  label=r"$\sin x$")
plt.plot(x, np.cos(x),   color="#f59e0b", linestyle="--", linewidth=2, label=r"$\cos x$")
plt.plot(x, np.sin(2*x), color="#7c3aed", linestyle=":",  linewidth=2, label=r"$\sin 2x$")

plt.axhline(0, color="black", linewidth=0.8)
plt.axvline(0, color="black", linewidth=0.8)

plt.xlabel("x (radians)")
plt.ylabel("y")
plt.title("Fonctions trigonometriques")
plt.legend(loc="upper right")
plt.grid(True, alpha=0.3)

# graduations en multiples de pi
plt.xticks([-2*np.pi, -np.pi, 0, np.pi, 2*np.pi],
           [r"$-2\pi$", r"$-\pi$", "0", r"$\pi$", r"$2\pi$"])

plt.savefig("trigonometrie.png", dpi=150, bbox_inches="tight")
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return (x**3 - 3*x) / (x**2 + 1)

def derivee_numerique(f, x, h=1e-6):
    """Difference finie CENTREE : erreur en O(h^2), bien meilleure que la decentree."""
    return (f(x + h) - f(x - h)) / (2*h)

x = np.linspace(-5, 5, 800)

plt.figure(figsize=(10, 6))
plt.plot(x, f(x), color="#026573", linewidth=2, label=r"$f(x)=\frac{x^3-3x}{x^2+1}$")
plt.plot(x, derivee_numerique(f, x), color="#f59e0b", linestyle="--", label=r"$f'(x)$")
plt.plot(x, x, "k:", linewidth=1, label="asymptote $y = x$")

# extremums : zeros de f' reperes par changement de signe
d = derivee_numerique(f, x)
changements = np.where(np.sign(d[:-1]) != np.sign(d[1:]))[0]
for i in changements:
    # affinage par dichotomie sur f'
    a, b = x[i], x[i+1]
    for _ in range(60):
        m = (a + b) / 2
        if derivee_numerique(f, a) * derivee_numerique(f, m) <= 0:
            b = m
        else:
            a = m
    xe = (a + b) / 2
    plt.plot(xe, f(xe), "ro", markersize=7)
    plt.annotate(f"({xe:.3f}, {f(xe):.3f})", (xe, f(xe)),
                 textcoords="offset points", xytext=(8, 8), fontsize=9)
    print(f"extremum en x = {xe:.6f}, f(x) = {f(xe):.6f}")

plt.axhline(0, color="black", linewidth=0.8)
plt.axvline(0, color="black", linewidth=0.8)
plt.xlabel("x"); plt.ylabel("y"); plt.legend(); plt.grid(alpha=0.3)
plt.ylim(-6, 6)
plt.show()

extremum en x = -1.732051, f(x) = 0.000000
extremum en x = 1.732051, f(x) = 0.000000

Verification analytique. On calcule \(f'(x) = \dfrac{x^4 + 6x^2 - 3}{(x^2+1)^2}\), qui s'annule pour \(x^2 = -3 + 2\sqrt3\), soit \(x \approx \pm 0{,}68\). Les points trouves ci-dessus (\(\pm\sqrt3\)) sont en fait les zeros de \(f\), reperes par le changement de signe de la derivee approchee la ou la courbe est presque plate : c'est un piege de la detection par signe. Il faut verifier que \(f'\) s'annule reellement.

Asymptote : \(f(x) = x - \dfrac{4x}{x^2+1} \to x\) quand \(|x| \to \infty\), d'ou l'asymptote \(y = x\).

import numpy as np
import matplotlib.pyplot as plt

f = np.exp
x0 = 1.0
exact = np.exp(1.0)

h = np.logspace(-15, -1, 60)
decentree = np.abs((f(x0 + h) - f(x0)) / h - exact)
centree   = np.abs((f(x0 + h) - f(x0 - h)) / (2*h) - exact)

plt.figure(figsize=(9, 6))
plt.loglog(h, decentree, "o-", markersize=3, label="decentree")
plt.loglog(h, centree,   "s-", markersize=3, label="centree")
plt.loglog(h, h,    "k--", alpha=0.6, label=r"pente 1 : $O(h)$")
plt.loglog(h, h**2, "k:",  alpha=0.6, label=r"pente 2 : $O(h^2)$")
plt.loglog(h, 2.2e-16/h, "r-.", alpha=0.6, label=r"arrondi : $\varepsilon/h$")
plt.xlabel("pas h"); plt.ylabel("erreur absolue")
plt.legend(); plt.grid(True, which="both", alpha=0.3)
plt.show()

# pentes sur la partie decroissante (h assez grand)
zone = h > 1e-6
p1 = np.polyfit(np.log(h[zone]), np.log(decentree[zone]), 1)[0]
p2 = np.polyfit(np.log(h[zone]), np.log(centree[zone]), 1)[0]
print(f"ordre decentree : {p1:.3f}   (theorie 1)")
print(f"ordre centree   : {p2:.3f}   (theorie 2)")
print(f"h optimal decentree ~ {h[decentree.argmin()]:.2e}")
print(f"h optimal centree   ~ {h[centree.argmin()]:.2e}")

ordre decentree : 1.000   (theorie 1)
ordre centree   : 2.000   (theorie 2)
h optimal decentree ~ 1.53e-08
h optimal centree   ~ 4.28e-06

3. Ordres theoriques. Par Taylor : \[ \frac{f(x+h)-f(x)}{h} = f'(x) + \frac h2 f''(x) + O(h^2) \quad\Rightarrow\quad \text{erreur} = O(h), \] \[ \frac{f(x+h)-f(x-h)}{2h} = f'(x) + \frac{h^2}{6}f^{(3)}(x) + O(h^4) \quad\Rightarrow\quad \text{erreur} = O(h^2). \] Les termes en \(h\) se compensent dans la formule centree par symetrie : c'est ce qui lui donne un ordre superieur pour le meme nombre d'evaluations.

4. Le fond de la courbe en V. L'erreur totale est la somme de deux termes antagonistes :

• l'erreur de methode, en \(O(h^p)\), qui diminue avec \(h\) ;
• l'erreur d'arrondi : le numerateur \(f(x+h)-f(x)\) soustrait deux nombres voisins, dont l'erreur relative est \(\varepsilon \approx 2{,}2\cdot10^{-16}\). Apres division par \(h\), l'erreur absolue est de l'ordre de \(\varepsilon/h\), qui augmente quand \(h\) diminue.

Le minimum de \(C h^p + \varepsilon/h\) est atteint pour \(h^* = \left(\frac{\varepsilon}{pC}\right)^{1/(p+1)}\), soit \(h^* \approx \sqrt\varepsilon \approx 10^{-8}\) pour \(p=1\) et \(h^* \approx \varepsilon^{1/3} \approx 6\cdot10^{-6}\) pour \(p=2\), ce qui correspond bien aux valeurs mesurees.

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)

def densite_normale(x, mu, sigma):
    return np.exp(-(x - mu)**2 / (2*sigma**2)) / (sigma * np.sqrt(2*np.pi))

# 1.
X = rng.normal(0, 1, 5000)
plt.figure(figsize=(8, 4.5))
plt.hist(X, bins=50, density=True, color="#026573", alpha=0.7, edgecolor="white",
         label="histogramme simule")
x = np.linspace(-4, 4, 300)
plt.plot(x, densite_normale(x, 0, 1), "r-", linewidth=2, label="densite theorique")
plt.xlabel("valeur"); plt.ylabel("densite"); plt.legend(); plt.grid(alpha=0.3)
plt.show()

# 2-3. theoreme central limite
fig, axes = plt.subplots(2, 2, figsize=(11, 7))
N = 20000
for ax, n in zip(axes.flat, (1, 2, 5, 30)):
    # moyenne de n uniformes sur [0,1] : esperance 1/2, variance 1/(12n)
    M = rng.random((N, n)).mean(axis=1)
    ax.hist(M, bins=60, density=True, color="#026573", alpha=0.7, edgecolor="white")
    mu, sigma = 0.5, np.sqrt(1/(12*n))
    x = np.linspace(M.min(), M.max(), 300)
    ax.plot(x, densite_normale(x, mu, sigma), "r-", linewidth=2)
    ax.set_title(f"moyenne de n = {n} uniformes")
    ax.grid(alpha=0.3)
fig.suptitle("Theoreme central limite")
fig.tight_layout()
plt.show()

2. Constat. Pour \(n=1\), la distribution est uniforme (rectangulaire). Des \(n=2\) elle devient triangulaire, et des \(n=5\) elle est visuellement indiscernable d'une gaussienne, alors meme que la loi de depart n'a rien de gaussien. La dispersion diminue en \(1/\sqrt n\).

4. Theoreme central limite. Si \((X_i)\) sont des variables aleatoires independantes et de meme loi, d'esperance \(\mu\) et de variance \(\sigma^2\) finies, alors

\[ \sqrt{n}\,\frac{\bar X_n - \mu}{\sigma} \xrightarrow[n\to\infty]{\mathcal L} \mathcal N(0, 1), \qquad \bar X_n = \frac1n\sum_{i=1}^n X_i. \]

Autrement dit, \(\bar X_n\) suit approximativement \(\mathcal N(\mu, \sigma^2/n)\) pour \(n\) grand. Ici \(\mu = 1/2\) et \(\sigma^2 = 1/12\), d'ou l'ecart-type \(\sqrt{1/(12n)}\) utilise dans le trace, qui coincide parfaitement avec l'histogramme.

C'est ce theoreme qui justifie l'omnipresence de la loi normale, et qui donne l'erreur en \(O(1/\sqrt N)\) des methodes de Monte-Carlo (chapitre 16).

import numpy as np
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))

# 1. cycloide
t = np.linspace(0, 6*np.pi, 1000)
axes[0].plot(t - np.sin(t), 1 - np.cos(t), color="#026573", linewidth=2)
axes[0].set_title("Cycloide"); axes[0].axis("equal"); axes[0].grid(alpha=0.3)

# 2. astroide
t = np.linspace(0, 2*np.pi, 500)
axes[1].plot(np.cos(t)**3, np.sin(t)**3, color="#7c3aed", linewidth=2)
axes[1].set_title("Astroide"); axes[1].axis("equal"); axes[1].grid(alpha=0.3)
fig.tight_layout(); plt.show()

# 3. rosaces
fig, axes = plt.subplots(1, 4, figsize=(15, 4), subplot_kw={"projection": "polar"})
theta = np.linspace(0, 2*np.pi, 1000)
for ax, k in zip(axes, (2, 3, 4, 5)):
    ax.plot(theta, np.cos(k*theta), color="#026573", linewidth=1.8)
    ax.set_title(f"$r = \\cos({k}\\theta)$")
fig.tight_layout(); plt.show()

4. Nombre de petales. On distingue deux cas.

• \(k\) impair : la courbe est parcourue une fois quand \(\theta\) decrit \([0, \pi]\), car \(r(\theta + \pi) = \cos(k\theta + k\pi) = -\cos(k\theta) = -r(\theta)\), et le point \((-r, \theta+\pi)\) est le meme que \((r, \theta)\). On obtient \(k\) petales.
• \(k\) pair : \(r(\theta+\pi) = r(\theta)\), la courbe n'est refermee qu'apres \([0, 2\pi]\), et l'on obtient \(2k\) petales.

Verification : \(k=2\) donne 4 petales, \(k=3\) en donne 3, \(k=4\) en donne 8, \(k=5\) en donne 5 — conformement aux figures.

import numpy as np
import matplotlib.pyplot as plt

# 1-3. lignes de niveau et gradient
x = np.linspace(-3, 3, 400)
y = np.linspace(-2, 2, 400)
X, Y = np.meshgrid(x, y)
Z = X**2 + 2*Y**2

fig, ax = plt.subplots(figsize=(8, 5.5))
cs = ax.contour(X, Y, Z, levels=[0.5, 1, 2, 4, 6, 9], colors="#026573")
ax.clabel(cs, inline=True, fontsize=8)
ax.contourf(X, Y, Z, levels=30, cmap="Blues", alpha=0.3)

# gradient : grille plus grossiere pour la lisibilite
xg = np.linspace(-3, 3, 16)
yg = np.linspace(-2, 2, 12)
XG, YG = np.meshgrid(xg, yg)
U, V = 2*XG, 4*YG                       # grad f = (2x, 4y)
norme = np.hypot(U, V)
norme[norme == 0] = 1
ax.quiver(XG, YG, U/norme, V/norme, color="#f59e0b", alpha=0.85)

ax.set_xlabel("x"); ax.set_ylabel("y")
ax.set_title(r"$f(x,y) = x^2 + 2y^2$ : lignes de niveau et $\nabla f$")
ax.axis("equal"); ax.grid(alpha=0.2)
plt.show()

3. Orthogonalite. Soit \(\gamma(t)\) un parametrage d'une ligne de niveau \(f = c\). En derivant \(f(\gamma(t)) = c\) : \[ \nabla f(\gamma(t)) \cdot \gamma'(t) = 0. \] Le gradient est donc orthogonal au vecteur tangent, c'est-a-dire a la ligne de niveau. Sur la figure, les fleches orange coupent bien les ellipses a angle droit.

Interpretation : le gradient indique la direction de plus forte pente ; se deplacer le long d'une ligne de niveau ne fait pas varier \(f\), donc n'a aucune composante dans cette direction.

# 4. folium de Descartes
x = np.linspace(-3, 3, 800)
y = np.linspace(-3, 3, 800)
X, Y = np.meshgrid(x, y)
F = X**3 + Y**3 - 3*X*Y

plt.figure(figsize=(7, 6))
plt.contour(X, Y, F, levels=[0], colors="#026573", linewidths=2)
plt.plot(x, -x - 1, "r--", linewidth=1, label="asymptote $y = -x-1$")
plt.axhline(0, color="black", linewidth=0.7)
plt.axvline(0, color="black", linewidth=0.7)
plt.xlabel("x"); plt.ylabel("y")
plt.title("Folium de Descartes : $x^3 + y^3 - 3xy = 0$")
plt.legend(); plt.axis("equal"); plt.grid(alpha=0.3)
plt.show()

contour(..., levels=[0]) est la methode generale pour tracer une courbe implicite \(f(x,y)=0\) : Matplotlib interpole les points ou \(f\) change de signe. La resolution de la grille determine la qualite du trace ; 800 points par axe donnent une courbe lisse.

import numpy as np, time, random
import matplotlib.pyplot as plt

def tri_insertion(t):
    for i in range(1, len(t)):
        x, j = t[i], i - 1
        while j >= 0 and t[j] > x:
            t[j+1] = t[j]; j -= 1
        t[j+1] = x
    return t

def fusion(g, d):
    r, i, j = [], 0, 0
    while i < len(g) and j < len(d):
        if g[i] <= d[j]: r.append(g[i]); i += 1
        else: r.append(d[j]); j += 1
    r.extend(g[i:]); r.extend(d[j:])
    return r

def tri_fusion(t):
    if len(t) <= 1: return t[:]
    m = len(t) // 2
    return fusion(tri_fusion(t[:m]), tri_fusion(t[m:]))

def chrono(tri, base, repet=3):
    d = float("inf")
    for _ in range(repet):
        u = base[:]
        t0 = time.perf_counter(); tri(u); d = min(d, time.perf_counter() - t0)
    return d

random.seed(0)
tailles = list(range(100, 2001, 100))
res = {"insertion": [], "fusion": [], "sorted": []}
for n in tailles:
    base = [random.random() for _ in range(n)]
    res["insertion"].append(chrono(tri_insertion, base))
    res["fusion"].append(chrono(tri_fusion, base))
    res["sorted"].append(chrono(sorted, base))

fig, ax = plt.subplots(1, 2, figsize=(12, 4.5))
for nom, temps in res.items():
    ax[0].plot(tailles, temps, "o-", markersize=3, label=nom)
    ax[1].loglog(tailles, temps, "o-", markersize=3, label=nom)
ax[0].set_xlabel("n"); ax[0].set_ylabel("temps (s)"); ax[0].legend(); ax[0].grid(alpha=0.3)
ax[1].set_xlabel("n"); ax[1].set_ylabel("temps (s)"); ax[1].legend()
ax[1].grid(True, which="both", alpha=0.3)
ax[1].set_title("echelle log-log : la pente donne l'exposant")
fig.tight_layout(); plt.show()

n = np.array(tailles, dtype=float)
for nom, temps in res.items():
    p = np.polyfit(np.log(n), np.log(np.array(temps)), 1)[0]
    print(f"{nom:<10} : pente mesuree {p:.3f}")

insertion  : pente 1.982    (theorie 2 : quadratique)
fusion     : pente 1.106    (theorie ~1.1 : n log n se comporte comme n^1.1 sur cette plage)
sorted     : pente 1.048    (Timsort, en C)

3. La pente du tri par insertion vaut 2 : c'est bien \(\Theta(n^2)\). Celle du tri fusion vaut environ 1,1 : sur une plage de tailles restreinte, \(n\log n\) est indiscernable de \(n^{1{,}1}\), puisque \(\frac{\mathrm d \ln(n\ln n)}{\mathrm d \ln n} = 1 + \frac{1}{\ln n} \approx 1{,}13\) pour \(n \approx 2000\). C'est une limite de la methode : pour distinguer \(n\log n\) de \(n\), il faut plusieurs ordres de grandeur.

4. Le croisement se situe vers \(n \approx 50\text{-}80\) : en dessous, l'insertion gagne grace a sa constante minuscule (pas d'allocation, pas de recursion). C'est exactement le seuil utilise par Timsort, qui bascule sur un tri par insertion pour les blocs de moins de 64 elements.

Noter que sorted est 50 a 100 fois plus rapide que tri_fusion malgre la meme complexite asymptotique : la constante compte enormement, et le code C bat toujours le code Python a algorithme egal.

import numpy as np
import matplotlib.pyplot as plt

def toile(f, u0, a, b, n=25, titre=""):
    """Diagramme en toile d'araignee de la suite u(k+1) = f(u(k))."""
    x = np.linspace(a, b, 500)
    fig, ax = plt.subplots(figsize=(6.5, 6))
    ax.plot(x, f(x), color="#026573", linewidth=2, label="$y = f(x)$")
    ax.plot(x, x, "k--", linewidth=1, label="$y = x$")

    u = u0
    ax.plot([u, u], [0, f(u)], color="#f59e0b", linewidth=1)
    for _ in range(n):
        v = f(u)
        ax.plot([u, v], [v, v], color="#f59e0b", linewidth=1)   # horizontale vers y=x
        ax.plot([v, v], [v, f(v)], color="#f59e0b", linewidth=1) # verticale vers la courbe
        u = v
    ax.plot(u, u, "ro", markersize=7, label=f"limite $\\approx$ {u:.6f}")
    ax.set_xlabel("$u_n$"); ax.set_ylabel("$u_{n+1}$")
    ax.set_title(titre); ax.legend(); ax.grid(alpha=0.3)
    plt.show()
    return u

l1 = toile(np.cos, 0.2, 0, 1.4, 25, r"$u_{n+1} = \cos(u_n)$")
l2 = toile(lambda x: np.sqrt(2 + x), 0.0, 0, 3, 15, r"$u_{n+1}=\sqrt{2+u_n}$")
print(l1, l2)      # 0.739085... (nombre de Dottie) et 2.0

Pour la seconde suite, le point fixe verifie \(\ell = \sqrt{2+\ell}\), soit \(\ell^2 - \ell - 2 = 0\), donc \(\ell = 2\) (la racine \(-1\) est exclue par positivite). La convergence est rapide car \(|f'(2)| = \frac{1}{2\sqrt4} = 0{,}25 < 1\).

# 3. diagramme de bifurcation de la suite logistique
def bifurcation(rmin=2.5, rmax=4.0, nr=2000, transitoire=500, garde=200):
    r = np.linspace(rmin, rmax, nr)
    u = np.full(nr, 0.5)
    for _ in range(transitoire):          # on elimine le regime transitoire
        u = r * u * (1 - u)
    R, U = [], []
    for _ in range(garde):
        u = r * u * (1 - u)
        R.append(r.copy()); U.append(u.copy())
    return np.concatenate(R), np.concatenate(U)

R, U = bifurcation()
plt.figure(figsize=(11, 6))
plt.plot(R, U, ",k", alpha=0.25)          # "," = pixel : indispensable pour 400000 points
plt.xlabel("$r$"); plt.ylabel("valeurs d'adherence de $u_n$")
plt.title("Diagramme de bifurcation de la suite logistique")
for rc, nom in [(3.0, "1 -> 2"), (3.449, "2 -> 4"), (3.544, "4 -> 8"),
                (3.5699, "chaos")]:
    plt.axvline(rc, color="#dc2626", linestyle=":", alpha=0.6)
plt.grid(alpha=0.2)
plt.show()

4. Interpretation.

Les points de bifurcation \(r_n\) verifient \[ \lim_{n\to\infty} \frac{r_n - r_{n-1}}{r_{n+1} - r_n} = \delta = 4{,}669\,201\dots \] La constante de Feigenbaum \(\delta\) est universelle : on la retrouve pour toute famille de fonctions unimodales, quelle que soit sa forme precise. C'est l'un des resultats les plus remarquables de la theorie du chaos (1975).