Module L12 · Partie G · Mathématiques de l’informatique

Calcul numérique : dériver, intégrer, simuler, optimiser.

Les équations de la physique n’ont presque jamais de solution exacte ; l’ordinateur les résout approximativement, et l’art est de contrôler l’erreur. Ce module construit les outils du simulateur et de l’optimiseur : différences finies, intégration numérique, Euler et Runge-Kutta pour les équations différentielles, Newton pour les racines, descente de gradient pour les minima. Tout ce que fait un moteur physique, un contrôleur ou un entraînement de réseau de neurones.

Durée : 3 séances · Prérequis : L10, dérivées et primitives (1re/Tle), séance 15. Objectifs : dérivation numérique et son erreur, intégration (rectangles, trapèzes, Simpson, Monte-Carlo), EDO (Euler explicite/implicite, RK4, stabilité, pas adaptatif), équations non linéaires (dichotomie, Newton), optimisation (gradient, momentum, Newton, contraintes simples), différentiation automatique.

Ce que vous saurez faire à la fin

Simuler un pendule, un ressort, une orbite, un robot avec RK4 et savoir pourquoi Euler explose.
Minimiser une fonction de 100 variables par descente de gradient et diagnostiquer la convergence.
Écrire une différentiation automatique minimale (l’idée de PyTorch).
Estimer l’ordre d’une méthode numérique par l’expérience.

Références : Numerical Recipes, Strang Computational Science and Engineering, cours MIT 18.085, programme de physique-informatique MP.

Fiche de cours · Définitions

Calcul numérique : définitions

Définition (nombre flottant IEEE 754). x = ±(1,m)₂ × 2^e avec une mantisse de 52 bits (double) : 53 bits de précision, ε_machine = 2⁻⁵² ≈ 2,2·10⁻¹⁶ (plus petit ε tel que 1 + ε ≠ 1). Tout réel x est représenté par fl(x) = x(1 + δ), |δ| ≤ ε/2. Valeurs spéciales : ±inf, nan, −0. Les entiers jusqu’à 2⁵³ sont exacts.

Définition (erreur, conditionnement, stabilité). Erreur absolue |x̂ − x|, relative |x̂ − x|/|x|. Le conditionnement d’un problème mesure l’amplification des erreurs sur les données par la solution exacte (propriété du problème). Un algorithme est stable si ses erreurs d’arrondi sont du même ordre que celles causées par une petite perturbation des données (propriété de l’algorithme).

Définition (ordre de convergence). Une suite x_k → x* converge à l’ordre p si |x_k+1 − x*| ≤ C|x_k − x*|^p. Linéaire (p = 1, C < 1) : les chiffres corrects croissent d’une constante par itération ; quadratique (p = 2) : ils doublent.

Définition (ordre d’une méthode numérique). Une méthode d’intégration ou d’EDO est d’ordre p si son erreur globale est O(h^p), h étant le pas. Diviser h par 2 divise l’erreur par 2^p.

Définition (problème raide, stabilité d’un schéma). Une EDO est raide si elle contient des échelles de temps très différentes. Un schéma est A-stable s’il reste borné sur y′ = λy, Re λ < 0, pour tout pas h : Euler implicite l’est, Euler explicite exige |1 + hλ| < 1.

Fiche de cours · Formules

Formules à connaître

Méthode	Formule	Erreur
Différence avant	f′(x) ≈ (f(x+h) − f(x))/h	O(h) + ε/h → h optimal ≈ √ε ≈ 10⁻⁸
Différence centrée	f′(x) ≈ (f(x+h) − f(x−h))/2h	O(h²) + ε/h → h ≈ ε^1/3 ≈ 6·10⁻⁶
Rectangles / point milieu	h·Σ f(x_i + h/2)	O(h²)
Trapèzes	h·(f₀/2 + f₁ + … + f_n−1 + f_n/2)	O(h²), erreur = −(b−a)h²f″(ξ)/12
Simpson	(h/3)·(f₀ + 4f₁ + 2f₂ + 4f₃ + … + f_n)	O(h⁴)
Euler explicite	y_n+1 = y_n + h·f(t_n, y_n)	O(h)
Euler implicite	y_n+1 = y_n + h·f(t_n+1, y_n+1)	O(h), A-stable
RK4	k₁ = f(t, y), k₂ = f(t+h/2, y+hk₁/2), k₃ = f(t+h/2, y+hk₂/2), k₄ = f(t+h, y+hk₃) ; y + h(k₁+2k₂+2k₃+k₄)/6	O(h⁴)
Verlet (x″ = a(x))	x_n+1 = 2x_n − x_n−1 + h²a(x_n)	O(h²), conserve l’énergie

Newton : x_k+1 = x_k − f(x_k)/f′(x_k) ; en dimension n : x_k+1 = x_k − J(x_k)⁻¹F(x_k)

Descente de gradient : x_k+1 = x_k − α∇f(x_k) ; converge pour α < 2/L (L = constante de Lipschitz du gradient = plus grande valeur propre de la hessienne)

Conditionnement de f en x : κ = |x f′(x)/f(x)| ; d’une matrice : κ(A) = ‖A‖‖A⁻¹‖ = σ₁/σ_n on perd log₁₀ κ chiffres significatifs

Fiche de cours · Théorèmes et démonstrations

Démonstrations à savoir refaire (1/2)

Théorème 1 (convergence quadratique de Newton). Si f est C², f(x*) = 0, f′(x*) ≠ 0, alors pour x₀ assez proche de x*, la suite de Newton vérifie |x_k+1 − x*| ≤ C|x_k − x*|².

Taylor à l’ordre 2 autour de x_k : 0 = f(x*) = f(x_k) + f′(x_k)(x* − x_k) + ½f″(ξ)(x* − x_k)². Divisons par f′(x_k) : 0 = f(x_k)/f′(x_k) + x* − x_k + f″(ξ)(x* − x_k)²/(2f′(x_k)). Or x_k+1 = x_k − f(x_k)/f′(x_k), donc x_k+1 − x* = f″(ξ)(x_k − x*)²/(2f′(x_k)). Près de x*, |f″/(2f′)| ≤ C. Donc l’erreur est au carré à chaque pas : le nombre de chiffres corrects double. Si f′(x*) = 0 (racine double), la convergence redevient linéaire.

Théorème 2 (erreur de la méthode des trapèzes). Pour f C² sur [a, b] et n sous-intervalles de largeur h, |∫f − T_h| ≤ (b − a)h² max|f″|/12.

Sur un sous-intervalle [x_i, x_i+1], l’interpolation linéaire p vérifie f(x) − p(x) = f″(ξ)(x − x_i)(x − x_i+1)/2 (reste de l’interpolation de Lagrange). Intégrer : ∫(x − x_i)(x − x_i+1)dx = −h³/6, donc l’erreur locale est −h³f″(ξ_i)/12. Sommer les n = (b − a)/h morceaux : erreur globale = −(b − a)h²f″(ξ)/12 (moyenne des f″ par le théorème des valeurs intermédiaires). L’ordre 2 apparaît parce que h³ × n = (b − a)h².

Théorème 3 (Euler explicite est d’ordre 1). Pour y′ = f(t, y) avec f lipschitzienne en y (constante L), l’erreur globale d’Euler vérifie |y(t_n) − y_n| ≤ (e^LT − 1)·h·max|y″|/(2L).

Erreur locale (Taylor) : y(t_n+1) = y(t_n) + hf(t_n, y(t_n)) + h²y″(ξ)/2, donc un pas exact commet O(h²). Propagation : e_n+1 = |y(t_n+1) − y_n+1| ≤ |y(t_n) − y_n| + h|f(t_n, y(t_n)) − f(t_n, y_n)| + h²M/2 ≤ (1 + hL)e_n + h²M/2. Récurrence : e_n ≤ (h²M/2)·Σ_k<n(1 + hL)^k = (hM/2L)((1 + hL)ⁿ − 1) ≤ (hM/2L)(e^LT − 1) avec nh = T. Erreur locale O(h²) sur T/h pas → globale O(h). Même raisonnement : RK4 a une erreur locale O(h⁵), globale O(h⁴).

Fiche de cours · Théorèmes et démonstrations

Démonstrations à savoir refaire (2/2)

Théorème 4 (stabilité d’Euler explicite sur y′ = λy). La suite y_n+1 = (1 + hλ)y_n reste bornée si et seulement si |1 + hλ| ≤ 1 ; pour λ < 0 réel, cela impose h ≤ 2/|λ|. Euler implicite y_n+1 = y_n/(1 − hλ) est borné pour tout h > 0.

y_n = (1 + hλ)ⁿy₀ ; une suite géométrique est bornée ssi la raison est de module ≤ 1. Pour λ réel négatif : −1 ≤ 1 + hλ ≤ 1 ⇔ h ≤ 2/|λ|. Implicite : raison 1/(1 − hλ), de module < 1 dès que Re λ < 0. Conséquence pratique : un système raide (λ = −10⁶, une constante de temps de 1 µs dans un circuit ou une raideur élevée) impose à Euler explicite h < 2·10⁻⁶ même si l’on ne s’intéresse qu’à la dynamique lente ; l’implicite permet un grand pas au prix de la résolution d’un système à chaque pas.

Théorème 5 (annulation catastrophique). Soustraire deux nombres proches perd des chiffres significatifs : si x et y sont connus à ε près relativement et |x − y| ≈ δ|x|, l’erreur relative sur x − y est ≈ ε/δ.

fl(x) = x(1 + δ₁), fl(y) = y(1 + δ₂) ; fl(x) − fl(y) = (x − y) + xδ₁ − yδ₂ ; l’erreur absolue ≈ |x|ε reste, mais le résultat vaut ≈ δ|x| : erreur relative ε/δ. Exemples : (1 − cos x)/x² pour x petit (utiliser 2sin²(x/2)/x²) ; les racines de ax² + bx + c par −b ± √Δ quand b² ≫ 4ac (calculer la racine « sûre » puis l’autre par x₁x₂ = c/a) ; la variance par E[X²] − E[X]² (préférer l’algorithme de Welford).

Théorème 6 (descente de gradient sur une fonction quadratique). Pour f(x) = ½xᵀAx − bᵀx avec A symétrique définie positive de valeurs propres 0 < μ ≤ … ≤ L, la descente x_k+1 = x_k − α∇f(x_k) converge pour 0 < α < 2/L, avec un taux max(|1 − αμ|, |1 − αL|) ; le meilleur α = 2/(μ + L) donne un taux (κ − 1)/(κ + 1), κ = L/μ.

∇f = Ax − b, x* = A⁻¹b. x_k+1 − x* = (I − αA)(x_k − x*). Dans la base propre de A, chaque composante est multipliée par 1 − αλ_i ; convergence ssi |1 − αλ_i| < 1 pour tous les λ_i, i.e. α < 2/L. Le taux est le pire facteur ; il est minimisé quand 1 − αμ = −(1 − αL), soit α = 2/(μ + L). Un mauvais conditionnement (κ grand) donne un taux proche de 1 : convergence lente en zigzag — motivation du momentum, d’Adam, et du préconditionnement.

Fiche de cours · Méthodes

Méthodes et pièges

Méthode — vérifier l’ordre d’une méthode. Calculer l’erreur pour h, h/2, h/4 sur un cas à solution connue ; le rapport des erreurs doit tendre vers 2^p ; tracer log(erreur) contre log(h) : la pente est p. Un ordre inférieur au prévu signale un bug ou une solution pas assez régulière.

Méthode — choisir un solveur d’EDO. Problème lisse et non raide : RK4 ou RK45 adaptatif (scipy.integrate.solve_ivp). Raide (constantes de temps très différentes, chimie, circuits, contact) : méthode implicite (BDF, Radau). Mécanique conservative sur longue durée (orbites, pendules) : intégrateur symplectique (Verlet, saute-mouton). Temps réel embarqué : pas fixe, Euler semi-implicite ou RK2 (coût borné).

Méthode — résoudre f(x) = 0. Encadrer la racine (signe change) → dichotomie (sûre, lente : 1 bit par itération) → Newton depuis l’encadrement (rapide, peut diverger) ; ou Brent (scipy.optimize.brentq) qui combine les deux. Toujours vérifier |f(x̂)| et le nombre d’itérations.

Pièges : comparer des flottants avec == ; additionner de petits nombres à un grand (sommer du plus petit au plus grand, ou math.fsum) ; 0,1 + 0,2 ≠ 0,3 ; un pas h trop petit en différences finies (l’arrondi domine) ; Newton sur une racine double ou avec f′ ≈ 0 ; inverser une matrice mal conditionnée ; tester une méthode seulement sur un cas où elle est exacte (polynôme de bas degré).

Fiche de cours · Exercices corrigés

Exercices corrigés

Exercice 1. Calculer √2 par Newton sur f(x) = x² − 2 depuis x₀ = 1, trois itérations, et vérifier le doublement des chiffres corrects.

Correction. x_k+1 = x_k − (x_k² − 2)/(2x_k) = (x_k + 2/x_k)/2 (méthode de Héron). x₁ = 1,5 (erreur 8,6·10⁻²) ; x₂ = 1,41667 (2,5·10⁻³) ; x₃ = 1,414216 (2,1·10⁻⁶) ; x₄ = 1,41421356237 (1,6·10⁻¹²). Les exposants −2, −3, −6, −12 doublent : convergence quadratique, conforme au Théorème 1 avec C = f″/(2f′) = 1/(2√2) ≈ 0,35.

Exercice 2. On intègre ∫₀¹ e^x dx = e − 1 par trapèzes avec n = 4, 8, 16. Donner les erreurs attendues et le nombre de points pour atteindre 10⁻⁸. Même question avec Simpson.

Correction. Trapèzes : erreur ≈ (b − a)h²f″/12 avec f″ ≈ e sur [0, 1], moyenne (e − 1) ≈ 1,72 : ≈ 0,143h². n = 4 (h = 0,25) : 9·10⁻³ ; n = 8 : 2,2·10⁻³ ; n = 16 : 5,6·10⁻⁴ (rapport 4 à chaque doublement : ordre 2). Pour 10⁻⁸ : h² ≈ 7·10⁻⁸, h ≈ 2,6·10⁻⁴, n ≈ 3 800 points. Simpson (ordre 4, erreur ≈ (b − a)h⁴f⁽⁴⁾/180 ≈ 0,0095h⁴) : n = 4 : 3,7·10⁻⁵ ; pour 10⁻⁸ : h ≈ 0,032, n ≈ 32 points. D’où l’intérêt des méthodes d’ordre élevé pour les fonctions lisses.

Exercice 3. Un robot suit ẋ = v cos θ, ẏ = v sin θ, θ̇ = ω avec v = 1 m/s, ω = 1 rad/s (cercle de rayon 1). On intègre sur 2π s avec Euler explicite à h = 0,1 puis 0,01. Prédire l’erreur de position finale et proposer mieux.

Correction. Solution exacte : retour au point de départ après 2π s. Euler explicite « sort » du cercle à chaque pas (il avance le long de la tangente) : le rayon croît d’un facteur √(1 + h²ω²) par pas, soit après N = 2π/h pas un rayon (1 + h²)^N/2 ≈ e^πh : h = 0,1 → rayon 1,37, erreur de position ≈ 0,4 m ; h = 0,01 → 1,032, erreur ≈ 0,03 m (ordre 1 : ÷10 quand h ÷10). Mieux : RK4 (erreur ≈ 10⁻⁵ à h = 0,1), ou — spécifique à ce modèle — l’intégration exacte de l’arc de cercle sur un pas : x += (v/ω)(sin(θ + ωh) − sin θ), y −= (v/ω)(cos(θ + ωh) − cos θ), θ += ωh, exacte pour v, ω constants sur le pas (c’est le modèle d’odométrie utilisé en L20).

01 / Dériver et intégrer

Différences finies : l’erreur a un ordre

import numpy as np
f, fprime = np.sin, np.cos          # on connaît la vraie dérivée pour mesurer l'erreur
x0 = 1.0
print(f"{'h':>10} {'avant':>12} {'centrée':>12}")
for h in [1e-1, 1e-2, 1e-3, 1e-4, 1e-6, 1e-8, 1e-10]:
    avant = (f(x0 + h) - f(x0)) / h                  # ordre 1 : erreur ∝ h
    centree = (f(x0 + h) - f(x0 - h)) / (2 * h)      # ordre 2 : erreur ∝ h²
    print(f"{h:10.0e} {abs(avant - fprime(x0)):12.2e} {abs(centree - fprime(x0)):12.2e}")
print("→ l'erreur DIMINUE avec h puis REMONTE : erreur d'arrondi ∝ ε/h. Optimum vers h ≈ √ε (avant) ou ε^(1/3) (centrée).")

# Gradient numérique d'une fonction de plusieurs variables : une différence centrée par composante
def gradient(f, x, h=1e-5):
    g = np.zeros_like(x)
    for i in range(len(x)):
        d = np.zeros_like(x); d[i] = h
        g[i] = (f(x + d) - f(x - d)) / (2 * h)
    return g
rosenbrock = lambda p: (1 - p[0])**2 + 100 * (p[1] - p[0]**2)**2
print(gradient(rosenbrock, np.array([0.5, 0.5])), "(exact : [-47, -50])")

Pourquoi centrée = ordre 2

Taylor : f(x ± h) = f(x) ± h f′(x) + h²/2 f″(x) ± h³/6 f‴(x) + … La différence avant laisse un reste en h/2 f″ ; la différence centrée annule les termes pairs et laisse h²/6 f‴. L’erreur d’arrondi vient de la soustraction de deux nombres proches (cancellation, module L10) : ε·|f|/h. Somme des deux : un optimum vers h ≈ 10⁻⁵ pour la centrée en double précision. Cette compétition troncature/arrondi est le fil rouge de tout le calcul numérique.

01 / Dériver et intégrer

Intégration numérique : rectangles, trapèzes, Simpson, Monte-Carlo

import numpy as np
f = lambda x: np.exp(-x**2)                     # pas de primitive élémentaire
exact = 0.7468241328124271                      # ∫₀¹ e^(−x²) dx

def rectangles(f, a, b, n):
    x = a + (np.arange(n) + 0.5) * (b - a) / n          # point milieu : ordre 2
    return (b - a) / n * f(x).sum()
def trapezes(f, a, b, n):
    x = np.linspace(a, b, n + 1); y = f(x)
    return (b - a) / n * (y[0] / 2 + y[1:-1].sum() + y[-1] / 2)
def simpson(f, a, b, n):                                  # n pair ; ordre 4
    x = np.linspace(a, b, n + 1); y = f(x)
    return (b - a) / (3 * n) * (y[0] + 4 * y[1:-1:2].sum() + 2 * y[2:-1:2].sum() + y[-1])
def monte_carlo(f, a, b, n, rng=np.random.default_rng(0)):
    return (b - a) * f(rng.uniform(a, b, n)).mean()      # erreur ∝ 1/√n, quelle que soit la dimension

print(f"{'n':>6} {'milieu':>10} {'trapèzes':>10} {'Simpson':>10} {'Monte-Carlo':>12}")
for n in (4, 16, 64, 256):
    print(f"{n:6d} {abs(rectangles(f,0,1,n)-exact):10.1e} {abs(trapezes(f,0,1,n)-exact):10.1e} {abs(simpson(f,0,1,n)-exact):10.1e} {abs(monte_carlo(f,0,1,n)-exact):12.1e}")
print("Quadrupler n : milieu/trapèzes ÷16 (ordre 2), Simpson ÷256 (ordre 4), Monte-Carlo ÷2 seulement.")

# Mais en dimension 10, Simpson demande n^10 points ; Monte-Carlo reste en 1/√n : il gagne
rng = np.random.default_rng(1); d = 10
X = rng.uniform(0, 1, (200_000, d))
print("volume de la boule unité en dim 10 :", (np.sum(X**2, axis=1) <= 1).mean() * 2**d, "(exact ≈ 2.55)")

Ordre p : l’erreur est en O(h^p). Pour vérifier l’ordre d’une méthode, doublez n et regardez de combien l’erreur est divisée : 2^p. C’est le test expérimental à faire systématiquement quand on implémente une méthode numérique.

02 / Équations différentielles

Euler et Runge-Kutta : simuler y′ = f(t, y)

import numpy as np
import matplotlib.pyplot as plt

def euler(f, y0, t):
    y = np.zeros((len(t), len(y0))); y[0] = y0
    for k in range(len(t) - 1):
        h = t[k + 1] - t[k]; y[k + 1] = y[k] + h * f(t[k], y[k])          # ordre 1
    return y

def rk4(f, y0, t):
    y = np.zeros((len(t), len(y0))); y[0] = y0
    for k in range(len(t) - 1):
        h = t[k + 1] - t[k]; tk, yk = t[k], y[k]
        k1 = f(tk, yk); k2 = f(tk + h / 2, yk + h / 2 * k1)
        k3 = f(tk + h / 2, yk + h / 2 * k2); k4 = f(tk + h, yk + h * k3)
        y[k + 1] = yk + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)                # ordre 4
    return y

# Oscillateur harmonique : x'' = −ω² x  ⇔  y = (x, v), y' = (v, −ω² x). Solution exacte : cos(ωt)
w = 2.0
f = lambda t, y: np.array([y[1], -w**2 * y[0]])
for h in (0.1, 0.05, 0.025):
    t = np.arange(0, 10 + h / 2, h)
    e_euler = abs(euler(f, [1, 0], t)[-1, 0] - np.cos(w * 10))
    e_rk4 = abs(rk4(f, [1, 0], t)[-1, 0] - np.cos(w * 10))
    print(f"h={h:<6} erreur Euler {e_euler:.2e}   RK4 {e_rk4:.2e}")

t = np.arange(0, 20, 0.05)
plt.figure(figsize=(8, 3))
plt.plot(t, euler(f, [1, 0], t)[:, 0], label="Euler h=0.05 (l'énergie augmente !)")
plt.plot(t, rk4(f, [1, 0], t)[:, 0], label="RK4"); plt.plot(t, np.cos(w * t), "--", label="exact")
plt.legend(); plt.show()

Pourquoi Euler « explose » et RK4 non

Euler explicite suit la tangente : sur un cercle (l’oscillateur), la tangente sort toujours vers l’extérieur, l’amplitude croît géométriquement. RK4 évalue la pente 4 fois par pas et combine (pondération de Simpson) : erreur locale en h⁵, globale en h⁴. Diviser h par 2 divise l’erreur d’Euler par 2, celle de RK4 par 16. Pour la mécanique (énergie à conserver sur des millions de pas — orbites, dynamique moléculaire), on préfère les schémas symplectiques (Verlet, saute-mouton) qui conservent une énergie approchée exactement. Sur un microcontrôleur, Euler avec un pas très petit reste courant par simplicité — à condition de connaître ses limites.

02 / Équations différentielles

Raideur, stabilité, pas adaptatif

import numpy as np
# Équation raide : y' = −100 (y − cos t). La solution suit cos t, mais Euler explicite exige h < 2/100.
f = lambda t, y: np.array([-100 * (y[0] - np.cos(t))])
for h in (0.03, 0.019, 0.005):
    t = np.arange(0, 3, h); y = euler(f, [0.0], t)
    print(f"Euler explicite h={h}: y(3) = {y[-1, 0]:.3g}  (exact ≈ {np.cos(3):.3f})")

# Euler implicite : y_{k+1} = y_k + h f(t_{k+1}, y_{k+1}) → inconditionnellement stable (ici résoluble à la main)
def euler_implicite(y0, t):
    y = [y0]
    for k in range(len(t) - 1):
        h = t[k + 1] - t[k]
        y.append((y[-1] + h * 100 * np.cos(t[k + 1])) / (1 + 100 * h))
    return np.array(y)
t = np.arange(0, 3, 0.03); print("Euler implicite h=0.03 :", euler_implicite(0.0, t)[-1].round(3))

# Pas adaptatif (idée de RK45 / scipy.integrate.solve_ivp) : comparer deux estimations, ajuster h
def rk4_pas(f, t, y, h):
    k1 = f(t, y); k2 = f(t + h/2, y + h/2*k1); k3 = f(t + h/2, y + h/2*k2); k4 = f(t + h, y + h*k3)
    return y + h/6*(k1 + 2*k2 + 2*k3 + k4)
def rk4_adaptatif(f, y0, t0, t1, tol=1e-6):
    t, y, h, pas = t0, np.array(y0, float), 0.1, 0
    while t < t1:
        h = min(h, t1 - t)
        y1 = rk4_pas(f, t, y, h); y2 = rk4_pas(f, t + h/2, rk4_pas(f, t, y, h/2), h/2)   # un pas h vs deux pas h/2
        err = np.abs(y2 - y1).max()
        if err < tol: t += h; y = y2; pas += 1
        h *= min(2, max(0.2, 0.9 * (tol / (err + 1e-16))**0.2))
    return y, pas
g = lambda t, y: np.array([y[1], -4 * y[0]])
y, pas = rk4_adaptatif(g, [1, 0], 0, 10); print(f"RK4 adaptatif : y(10) = {y[0]:.6f} (exact {np.cos(20):.6f}) en {pas} pas")

Une équation est raide quand elle contient des échelles de temps très différentes (une réaction chimique rapide et une lente, un moteur électrique et sa charge mécanique). Les méthodes explicites sont alors contraintes par l’échelle rapide même si on ne s’y intéresse pas ; les implicites résolvent un système à chaque pas mais acceptent de grands h. Sur PC : scipy.integrate.solve_ivp(f, (0, 10), y0, method="RK45") ou "LSODA" pour le raide.

02 / Équations différentielles

Application : le pendule inversé (le robot qui tient debout)

import numpy as np
import matplotlib.pyplot as plt
# Pendule inversé sur chariot, linéarisé près de la verticale : θ'' = (g/l) θ + u/(m l)
g, l = 9.81, 0.5
def pendule(u):                                    # u : fonction de commande (force) de l'état
    return lambda t, y: np.array([y[1], g / l * y[0] + u(y) / (0.2 * l)])
t = np.arange(0, 3, 0.005)

sans = rk4(pendule(lambda y: 0.0), [0.05, 0], t)                               # rien : il tombe
Kp, Kd = 20.0, 4.0
pd = rk4(pendule(lambda y: -Kp * y[0] - Kd * y[1]), [0.05, 0], t)              # régulateur PD (séance 29)
plt.figure(figsize=(8, 3))
plt.plot(t, sans[:, 0], label="sans commande"); plt.plot(t, pd[:, 0], label="PD : Kp=20, Kd=4")
plt.ylim(-0.3, 0.3); plt.ylabel("angle (rad)"); plt.legend(); plt.show()

# Les valeurs propres du système bouclé disent si c'est stable : y' = A y avec A = [[0, 1], [g/l − Kp/(m l), −Kd/(m l)]]
A = np.array([[0, 1], [g / l - Kp / (0.2 * l), -Kd / (0.2 * l)]])
print("valeurs propres bouclées :", np.linalg.eigvals(A).round(2), "→ parties réelles négatives = stable")
A0 = np.array([[0, 1], [g / l, 0]]); print("sans commande :", np.linalg.eigvals(A0).round(2), "→ une racine positive = instable")

Voilà le pont entre la simulation numérique (RK4), l’algèbre linéaire (valeurs propres du module L10) et le contrôle (module L21) : un système est stable si ses valeurs propres ont une partie réelle négative, et le rôle du régulateur est de les déplacer. Le Segway, la fusée qui atterrit, le robot bipède : tous des pendules inversés.

03 / Résoudre et optimiser

Racines : dichotomie sûre, Newton rapide

La dichotomie est linéaire (1 bit par itération) mais garantie ; Newton est quadratique mais demande la dérivée et un bon départ. En pratique : Brent (dichotomie + sécante) sur PC (scipy.optimize.brentq), Newton-Raphson dans les moteurs physiques et la cinématique inverse (module L10). En dimension n, Newton résout J·Δx = −F à chaque pas : un système linéaire.

03 / Résoudre et optimiser

Descente de gradient : le moteur de l’apprentissage

import numpy as np
import matplotlib.pyplot as plt
f = lambda p: (1 - p[0])**2 + 100 * (p[1] - p[0]**2)**2          # Rosenbrock : vallée courbe, minimum en (1, 1)
grad = lambda p: np.array([-2 * (1 - p[0]) - 400 * p[0] * (p[1] - p[0]**2), 200 * (p[1] - p[0]**2)])

def descente(grad, x0, lr, iters, momentum=0.0):
    x, v, traj = np.array(x0, float), np.zeros(2), []
    for _ in range(iters):
        v = momentum * v - lr * grad(x); x = x + v; traj.append(x.copy())
    return np.array(traj)

fig, ax = plt.subplots(figsize=(7, 4))
X, Y = np.meshgrid(np.linspace(-1.5, 1.5, 200), np.linspace(-0.5, 1.5, 200))
ax.contour(X, Y, np.log(1 + f([X, Y])), levels=25, linewidths=0.5)
for lr, mom, nom in [(0.001, 0, "gradient lr=1e-3"), (0.001, 0.9, "momentum 0.9"), (0.003, 0, "lr=3e-3 (instable)")]:
    tr = descente(grad, [-1, 1], lr, 2000, mom)
    ax.plot(tr[:, 0], tr[:, 1], ".-", ms=2, lw=0.7, label=f"{nom} → f={f(tr[-1]):.2g}")
ax.plot(1, 1, "w*", ms=12); ax.legend(fontsize=8); plt.show()

# Newton en 2D utilise la courbure (hessienne) : quelques itérations suffisent
def hess(p):
    x, y = p
    return np.array([[2 - 400 * y + 1200 * x**2, -400 * x], [-400 * x, 200]])
x = np.array([-1.0, 1.0])
for k in range(8):
    x = x - np.linalg.solve(hess(x), grad(x))
print("Newton :", x.round(6), "en 8 itérations")

Le pas d’apprentissage et ses variantes

Trop petit : lent. Trop grand : oscille ou diverge (la borne de stabilité est 2/L, L = plus grande valeur propre de la hessienne). Le momentum accumule la vitesse et traverse les vallées ; Adam (module L14) adapte le pas par coordonnée ; le gradient stochastique (SGD) utilise un sous-échantillon des données à chaque pas : bruité mais 1000 fois moins cher, et le bruit aide à sortir des minima locaux. Newton est imbattable en petite dimension mais la hessienne coûte O(n²) mémoire et O(n³) par pas : impossible pour un réseau à 10⁹ paramètres. Entre les deux : L-BFGS, gradient conjugué.

03 / Résoudre et optimiser

Différentiation automatique : l’idée de PyTorch en 40 lignes

import math

class Val:
    """Un nombre qui se souvient de comment il a été calculé, pour propager le gradient en arrière."""
    def __init__(self, v, parents=(), op=""):
        self.v, self.grad, self._parents, self._op, self._backward = v, 0.0, parents, op, lambda: None
    def __add__(self, o):
        o = o if isinstance(o, Val) else Val(o); r = Val(self.v + o.v, (self, o), "+")
        def back(): self.grad += r.grad; o.grad += r.grad             # d(a+b)/da = 1
        r._backward = back; return r
    def __mul__(self, o):
        o = o if isinstance(o, Val) else Val(o); r = Val(self.v * o.v, (self, o), "*")
        def back(): self.grad += o.v * r.grad; o.grad += self.v * r.grad   # d(ab)/da = b
        r._backward = back; return r
    def __pow__(self, n):
        r = Val(self.v ** n, (self,), f"**{n}")
        def back(): self.grad += n * self.v ** (n - 1) * r.grad
        r._backward = back; return r
    def tanh(self):
        t = math.tanh(self.v); r = Val(t, (self,), "tanh")
        def back(): self.grad += (1 - t * t) * r.grad
        r._backward = back; return r
    __radd__ = __add__; __rmul__ = __mul__
    def __neg__(self): return self * -1
    def __sub__(self, o): return self + (-o)
    def __rsub__(self, o): return Val(o) + (-self)
    def backward(self):
        ordre, vus = [], set()                    # tri topologique du graphe de calcul (module L05)
        def visiter(n):
            if n not in vus:
                vus.add(n)
                for p in n._parents: visiter(p)
                ordre.append(n)
        visiter(self); self.grad = 1.0
        for n in reversed(ordre): n._backward()   # règle de la chaîne, en arrière

x, y = Val(-1.0), Val(1.0)
f = (1 - x)**2 + 100 * (y - x**2)**2              # Rosenbrock, construit comme un graphe
f.backward()
print("f =", f.v, " ∂f/∂x =", x.grad, " ∂f/∂y =", y.grad, "(exact : -404, 0)")

# Un neurone : tanh(w1 x1 + w2 x2 + b), et le gradient par rapport aux poids
w1, w2, b = Val(0.5), Val(-0.3), Val(0.1)
sortie = (w1 * 2.0 + w2 * 1.0 + b).tanh()
sortie.backward()
print("sortie", round(sortie.v, 4), "gradients :", round(w1.grad, 4), round(w2.grad, 4), round(b.grad, 4))

Pourquoi c’est révolutionnaire

Les différences finies coûtent n évaluations pour n paramètres (et sont imprécises) ; la différentiation symbolique (module L02) explose en taille. La différentiation automatique en mode inverse (backpropagation) calcule le gradient complet en un passage arrière, au coût d’environ 2-3 évaluations de la fonction, quel que soit n. C’est ce qui rend possible l’entraînement de réseaux à des milliards de paramètres. PyTorch et JAX font exactement ceci, sur des tenseurs et sur GPU. Cette classe (inspirée de micrograd, Andrej Karpathy) est le point de départ du module L14.

Cours

Cours 1 — Consistance, stabilité, convergence : le théorème qui organise tout

Pour une méthode numérique d’EDO (ou d’intégration, de résolution) on distingue trois propriétés :

Consistance (ordre p) : l’erreur commise en un pas, à partir de la solution exacte, est O(h^p+1). Se vérifie par développement de Taylor. Euler : p = 1 ; RK4 : p = 4.
Stabilité : les erreurs (arrondi, troncature) ne sont pas amplifiées sans borne au fil des pas. Pour y′ = λy (λ < 0), Euler explicite exige |1 + hλ| < 1, soit h < 2/|λ| : c’est la région de stabilité. Euler implicite : stable pour tout h.
Convergence : l’erreur globale à t fixé tend vers 0 quand h → 0, en O(h^p).

Théorème de Lax (équivalence). Pour un problème linéaire bien posé, une méthode consistante converge si et seulement si elle est stable. Autrement dit : on ne peut pas juger une méthode par son ordre seul ; il faut la région de stabilité. C’est pourquoi une méthode d’ordre 4 explicite peut exploser là où Euler implicite d’ordre 1 tient.

Cours

Cours 2 — Exemple travaillé : une masse-ressort avec frottement, quatre schémas, mêmes données

import numpy as np
import matplotlib.pyplot as plt
# m x'' = −k x − c x'  ;  y = (x, v). Solution exacte connue : oscillation amortie
m, k, c = 1.0, 40.0, 0.4
f = lambda t, y: np.array([y[1], (-k * y[0] - c * y[1]) / m])
w0 = np.sqrt(k / m); zeta = c / (2 * np.sqrt(k * m)); wd = w0 * np.sqrt(1 - zeta**2)
exact = lambda t: np.exp(-zeta * w0 * t) * (np.cos(wd * t) + zeta * w0 / wd * np.sin(wd * t))     # x(0) = 1, v(0) = 0

def euler_expl(f, y0, t):
    y = [np.array(y0, float)]
    for a, b in zip(t, t[1:]): y.append(y[-1] + (b - a) * f(a, y[-1]))
    return np.array(y)
def euler_semi_impl(f, y0, t):                       # Euler symplectique : v d'abord, puis x avec le NOUVEAU v
    y = [np.array(y0, float)]
    for a, b in zip(t, t[1:]):
        h = b - a; x, v = y[-1]; v = v + h * (-k * x - c * v) / m; x = x + h * v; y.append(np.array([x, v]))
    return np.array(y)
def rk4_(f, y0, t):
    y = [np.array(y0, float)]
    for a, b in zip(t, t[1:]):
        h = b - a; yk = y[-1]; k1 = f(a, yk); k2 = f(a + h/2, yk + h/2*k1); k3 = f(a + h/2, yk + h/2*k2); k4 = f(b, yk + h*k3)
        y.append(yk + h/6 * (k1 + 2*k2 + 2*k3 + k4))
    return np.array(y)
h = 0.05; t = np.arange(0, 10 + h/2, h)
print(f"h = {h}, ω0 = {w0:.2f} → hω0 = {h * w0:.2f} (Euler explicite limite : 2 pour l'oscillateur... voir)")
for nom, m_ in [("Euler explicite", euler_expl), ("Euler semi-implicite", euler_semi_impl), ("RK4", rk4_)]:
    y = m_(f, [1, 0], t); err = np.abs(y[:, 0] - exact(t)).max()
    print(f"{nom:22s} erreur max {err:.2e}"); plt.plot(t, y[:, 0], label=nom)
plt.plot(t, exact(t), "w--", lw=1, label="exact"); plt.ylim(-1.5, 1.5); plt.legend(fontsize=8); plt.show()

Lecture. Euler explicite gagne de l’énergie à chaque pas (amplification |1 + z| > 1 pour z imaginaire pur) ; le semi-implicite (une ligne de différence : utiliser le nouveau v pour x) est symplectique et reste borné ; RK4 est précis mais dériverait sur des millions de pas. Sur un microcontrôleur qui simule ou filtre à 1 kHz, le semi-implicite est le meilleur rapport coût/robustesse — c’est ce qu’utilisent les moteurs physiques des jeux.

Cours

Cours 3 — Optimisation : conditions d’optimalité, convexité, contraintes

Notion	Énoncé	Usage
Condition du premier ordre	En un minimum intérieur, ∇f(x*) = 0	Newton résout ∇f = 0 ; le gradient s’arrête quand ‖∇f‖ < ε
Second ordre	Hessienne semi-définie positive en x*	Distingue minimum, maximum, point-selle (fréquent en grande dimension)
Convexité	f(λx + (1−λ)y) ≤ λf(x) + (1−λ)f(y) ⇔ hessienne ⪰ 0	Tout minimum local est global ; le gradient converge. Moindres carrés, régression logistique, LQR : convexes. Réseaux de neurones : non.
Lipschitz du gradient (L)	‖∇f(x) − ∇f(y)‖ ≤ L‖x − y‖	Pas de gradient sûr : lr < 2/L ; L = plus grande valeur propre de la hessienne
Forte convexité (μ)	hessienne ⪰ μI	Convergence linéaire au taux (1 − μ/L) : le conditionnement L/μ gouverne la vitesse
Contraintes d’égalité	Lagrangien L = f + Σ λ_i g_i ; ∇L = 0	Multiplicateurs de Lagrange
Contraintes d’inégalité	Conditions KKT ; projection ; barrières	MPC (L21), SVM, allocation de ressources

Avant d’optimiser, demandez-vous : la fonction est-elle convexe (alors le gradient suffit et le résultat est unique) ? Mal conditionnée (alors normaliser, préconditionner, ou Newton) ? Contrainte (alors projeter, pénaliser, ou un solveur QP) ? Ces trois questions couvrent 90 % des situations en IA et en robotique.

TP guidé

TP — Un simulateur de vol de drone 2D et son autopilote optimisé (sur PC, 3 h)

Modèle. Drone plan : état (x, z, θ, ẋ, ż, θ̇), commandes (poussée totale T, couple τ). Équations : m ẍ = −T sin θ, m z̈ = T cos θ − m g, I θ̈ = τ. Fichier drone.py avec f(t, y, u).
Intégrateurs. Euler, semi-implicite, RK4, et scipy.integrate.solve_ivp (RK45, pas adaptatif). Sans commande (T = mg, τ = 0) : le drone doit flotter ; avec τ constant : rotation ; vérifiez la conservation de l’énergie en chute libre. Mesurez l’ordre de chaque intégrateur sur un cas à solution exacte (chute libre).
Contrôle. Autopilote en cascade : PD sur la position → angle désiré → PD sur l’angle → couple ; PD sur l’altitude → poussée. Réglez à la main pour aller de (0, 0) à (5, 3) sans dépassement > 10 %.
Optimisation des gains. Coût = ∫ (erreur de position² + 0,01 commande²) dt sur une simulation de 8 s. Minimisez par scipy.optimize.minimize (Nelder-Mead, puis L-BFGS avec gradient par différences finies). Comparez au réglage manuel ; tracez la surface de coût pour deux gains.
Robustesse. Perturbez la masse de ±20 % et ajoutez un vent constant : les gains optimisés tiennent-ils ? Ajoutez un terme intégral et recommencez.
Livrable. Dépôt avec drone.py, integrateurs.py (+ test d’ordre), autopilote.py, optimiser.py, animation (matplotlib FuncAnimation ou GIF) et RESULTATS.md (ordres mesurés, gains, coûts, robustesse).

Exercices

Exercices auto-corrigés — dériver, intégrer, résoudre

Exercice 1 — Ordre d’une formule de dérivation

Écrivez derivee_ordre4(f, x, h) = (−f(x+2h) + 8f(x+h) − 8f(x−h) + f(x−2h)) / (12h) et ordre_experimental(methode, f, fprime, x) qui estime l’ordre par la pente log-log de l’erreur pour h = 10⁻¹, 10⁻², 10⁻³. Vérifiez : centrée → ≈ 2, ordre4 → ≈ 4.

Correction

def derivee_ordre4(f, x, h): return (-f(x + 2*h) + 8*f(x + h) - 8*f(x - h) + f(x - 2*h)) / (12*h)
def ordre_experimental(m, f, fp, x):
    hs = np.array([0.2, 0.1, 0.05, 0.025]); err = np.array([abs(m(f, x, h) - fp(x)) for h in hs])
    return np.polyfit(np.log(hs), np.log(err), 1)[0]

Exercice 2 — Intégration adaptative par Simpson récursif

simpson_adaptatif(f, a, b, tol) : calcule Simpson sur [a, b] et sur les deux moitiés ; si l’écart < 15·tol, renvoie la version raffinée, sinon récurse sur chaque moitié avec tol/2. Comptez les évaluations de f et comparez à Simpson à pas fixe pour la même précision sur ∫₀¹ √x dx (dérivée infinie en 0).

Correction

def simpson_adaptatif(f, a, b, tol=1e-8, fa=None, fb=None, fm=None):
    m = (a + b) / 2
    fa = f(a) if fa is None else fa; fb = f(b) if fb is None else fb; fm = f(m) if fm is None else fm
    S = (b - a) / 6 * (fa + 4 * fm + fb)
    fl, fr = f((a + m) / 2), f((m + b) / 2)
    Sl = (m - a) / 6 * (fa + 4 * fl + fm); Sr = (b - m) / 6 * (fm + 4 * fr + fb)
    if abs(Sl + Sr - S) < 15 * tol: return Sl + Sr + (Sl + Sr - S) / 15
    return simpson_adaptatif(f, a, m, tol / 2, fa, fm, fl) + simpson_adaptatif(f, m, b, tol / 2, fm, fb, fr)

Le facteur 15 vient de l’extrapolation de Richardson (l’erreur de Simpson est divisée par 16 quand h est divisé par 2). L’adaptatif concentre les points près de 0 où √x varie vite.

Exercices

Exercices auto-corrigés — EDO et optimisation

Exercice 3 — Verlet et énergie

Implémentez verlet(acc, x0, v0, h, n) (saute-mouton) pour un oscillateur x″ = −x et vérifiez que l’énergie ½v² + ½x² reste dans une bande étroite sur 10⁵ pas, alors que celle d’Euler explicite dérive.

Correction

def verlet(acc, x0, v0, h, n):
    x, v = np.empty(n + 1), np.empty(n + 1); x[0], v[0] = x0, v0; a = acc(x0)
    for k in range(n):
        vh = v[k] + h / 2 * a; x[k + 1] = x[k] + h * vh; a = acc(x[k + 1]); v[k + 1] = vh + h / 2 * a
    return x, v

Exercice 4 — Descente de gradient avec recherche linéaire

Implémentez gradient_armijo(f, grad, x0, iters) : à chaque pas, partir de α = 1 et diviser par 2 tant que f(x − α∇f) > f(x) − 10⁻⁴ α‖∇f‖² (condition d’Armijo). Sur Rosenbrock, elle doit approcher (1, 1) sans réglage de pas, là où un pas fixe de 0,01 diverge ou stagne.

Correction

def gradient_armijo(f, grad, x0, iters=5000):
    x = x0.copy()
    for _ in range(iters):
        g = grad(x); fx = f(x); alpha = 1.0
        while f(x - alpha * g) > fx - 1e-4 * alpha * g @ g: alpha /= 2
        x = x - alpha * g
    return x

04 / Défis

Défi ★ — Ordre expérimental et orbite

Consigne

1) Mesurez expérimentalement l’ordre d’Euler, du point milieu (RK2) et de RK4 sur y′ = −y (solution e^−t) : tracez l’erreur en fonction de h en échelle log-log et lisez la pente. 2) Simulez une orbite (gravitation : a = −GM r/|r|³) avec Euler puis RK4 pendant 20 périodes ; tracez l’énergie mécanique au cours du temps pour chaque méthode.

Correction (extrait)

def rk2(f, y0, t):
    y = np.zeros((len(t), len(y0))); y[0] = y0
    for k in range(len(t) - 1):
        h = t[k+1] - t[k]; k1 = f(t[k], y[k]); y[k+1] = y[k] + h * f(t[k] + h/2, y[k] + h/2 * k1)
    return y
f = lambda t, y: -y
hs = np.array([0.4, 0.2, 0.1, 0.05, 0.025])
for nom, m in [("Euler", euler), ("RK2", rk2), ("RK4", rk4)]:
    err = [abs(m(f, [1.0], np.arange(0, 2 + h/2, h))[-1, 0] - np.exp(-2)) for h in hs]
    pente = np.polyfit(np.log(hs), np.log(err), 1)[0]
    print(f"{nom}: ordre ≈ {pente:.2f}")

GM = 1.0
grav = lambda t, y: np.array([y[2], y[3], -GM * y[0] / np.hypot(y[0], y[1])**3, -GM * y[1] / np.hypot(y[0], y[1])**3])
t = np.arange(0, 20 * 2 * np.pi, 0.05)
for nom, m in [("Euler", euler), ("RK4", rk4)]:
    y = m(grav, [1, 0, 0, 1], t)
    E = 0.5 * (y[:, 2]**2 + y[:, 3]**2) - GM / np.hypot(y[:, 0], y[:, 1])
    print(nom, "dérive d'énergie :", f"{E[-1] - E[0]:+.3e}")

04 / Défis

Défi ★★ — Régression non linéaire par gradient et par Gauss-Newton

Consigne

Des mesures de décharge d’un condensateur : V(t) = V₀ e^−t/τ + bruit. 1) Estimez (V₀, τ) en minimisant la somme des carrés des résidus par descente de gradient (gradient par différences finies puis par votre classe Val). 2) Implémentez Gauss-Newton : linéariser les résidus (jacobienne J), résoudre JᵀJ Δ = −Jᵀ r à chaque itération. Comparez le nombre d’itérations. 3) Pourquoi le gradient est-il si lent ici ? (Regardez les échelles de V₀ et τ : conditionnement.)

Correction (extrait)

def residus(p): return p[0] * np.exp(-t / p[1]) - V
def jac(p):
    e = np.exp(-t / p[1])
    return np.column_stack([e, p[0] * t / p[1]**2 * e])
p = np.array([1.0, 0.5])
for k in range(10):
    J, r = jac(p), residus(p)
    p = p + np.linalg.solve(J.T @ J, -J.T @ r)
print("Gauss-Newton :", p.round(4), "en 10 itérations")
cout = lambda p: 0.5 * (residus(p)**2).sum()
p = np.array([1.0, 0.5])
for k in range(20000):
    p -= 0.002 * gradient(cout, p)
print("gradient :", p.round(4), "après 20 000 itérations")

Gauss-Newton exploite la structure « somme de carrés » (c’est Newton avec la hessienne approchée JᵀJ). Levenberg-Marquardt ajoute un amortissement λI pour la robustesse : c’est l’algorithme de scipy.optimize.curve_fit et de la calibration de caméras et de robots.

04 / Défis

Défi ★★★ — Esprit prépa : un moteur physique symplectique et son contrôleur

Consigne

1) Implémentez le schéma de Verlet (saute-mouton) pour un système de N masses reliées par des ressorts (chaîne) et montrez que l’énergie oscille sans dériver sur 10⁶ pas, contrairement à RK4 qui dérive lentement. 2) Prouvez que Verlet est d’ordre 2 et réversible en temps (en inversant les vitesses, on revient exactement au départ — vérifiez numériquement à 10⁻¹⁰ près). 3) Ajoutez au pendule inversé non linéaire (sin θ, pas θ) un régulateur trouvé par optimisation : minimisez par descente de gradient sur (Kp, Kd) une fonction de coût ∫(θ² + 0,01 u²) dt simulée ; comparez au PD choisi à la main. C’est un premier pas vers le contrôle optimal (module L21).

Correction (extrait : Verlet et réversibilité)

def verlet(acc, x0, v0, h, n):
    x, v = np.array(x0, float), np.array(v0, float)
    a = acc(x)
    for _ in range(n):
        v_demi = v + h / 2 * a
        x = x + h * v_demi
        a = acc(x)
        v = v_demi + h / 2 * a
    return x, v

k = 10.0
def acc_chaine(x):                 # N masses, ressorts entre voisines, extrémités fixes en 0
    xg = np.concatenate([[0], x, [0]])
    return k * (xg[:-2] - 2 * xg[1:-1] + xg[2:])
x0 = np.zeros(5); x0[2] = 0.3; v0 = np.zeros(5)
x, v = verlet(acc_chaine, x0, v0, 0.01, 100_000)
xr, vr = verlet(acc_chaine, x, -v, 0.01, 100_000)           # inverser le temps
print("retour au départ :", np.abs(xr - x0).max())         # ≈ 1e-11 : réversible
def energie(x, v):
    xg = np.concatenate([[0], x, [0]])
    return 0.5 * (v**2).sum() + 0.5 * k * (np.diff(xg)**2).sum()
print("énergie initiale", energie(x0, v0).round(6), "finale", energie(x, v).round(6))

Verlet est symplectique : il conserve exactement une énergie « voisine » de la vraie, donc pas de dérive séculaire. C’est le schéma de toute la dynamique moléculaire, des simulations astronomiques et de nombreux moteurs de jeux (Box2D, position-based dynamics). Pour la question 3, la descente de gradient sur les gains est une forme rudimentaire d’apprentissage par renforcement ; avec un modèle linéarisé, la solution exacte est le régulateur LQR (équation de Riccati).

05 / Vérification

Avec RK4 (ordre 4), diviser le pas par 2 divise l’erreur par environ :

Deux questions supplémentaires

1. Pourquoi ne pas prendre h = 10⁻¹⁵ pour une dérivée numérique ? L’erreur d’arrondi ε/h explose : la soustraction de deux nombres presque égaux ne laisse aucun chiffre significatif.

2. Quel avantage de la rétropropagation sur les différences finies pour 10⁶ paramètres ? Un seul passage arrière au lieu de 10⁶ évaluations.

Référence

Les mots à retenir

Mot	Définition
Ordre d’une méthode	p tel que l’erreur soit en O(h^p).
Troncature / arrondi	Erreur du schéma / erreur des flottants ; leur somme a un optimum.
Simpson	Intégration par paraboles, ordre 4.
Euler explicite	y + h f(t, y) ; ordre 1 ; stable seulement pour h petit.
RK4	4 pentes combinées ; ordre 4.
Raideur	Échelles de temps très différentes ; méthodes implicites.
Symplectique	Schéma qui conserve la structure hamiltonienne (Verlet).
Newton	x − f/f′ ; convergence quadratique.
Descente de gradient	x − lr·∇f ; momentum, SGD, Adam.
Gauss-Newton / LM	Moindres carrés non linéaires par linéarisation.
Différentiation automatique	Gradient exact par propagation arrière sur le graphe de calcul.

Pour continuer

Fin de la partie G : vous avez les maths de l’ingénieur

Partie H : l’intelligence artificielle. Module suivant : l’apprentissage supervisé from scratch — régression, classification, arbres, validation — avant les réseaux de neurones.

À faire chez soi

Refaire la séance 15 (simulation physique) avec RK4 et Verlet, et comparer.
Étendre la classe Val avec exp, log, la division, et vérifier chaque gradient par différences finies.
Lire le chapitre « Ordinary Differential Equations » de Numerical Recipes.

← L11 Sommaire L13 : Apprentissage supervisé →