Applications : algorithmes numériques de base

Cours · Python (école d'ingénieur), chapitre 10 · école d'ingénieur, informatique · le début du chapitre, la suite sur Adloun

10.1 Introduction et Motivation

Maintenant que nous maîtrisons les bases de Python ainsi que les structures conditionnelles, itératives, et les bibliothèques NumPy et Matplotlib, nous pouvons aborder nos premiers algorithmes de calcul numérique. Contrairement aux calculs sur le type int (précision arbitraire en Python), les calculs sur le type float comportent de nombreuses approximations dues au codage binaire. Ce chapitre explore comment concevoir des algorithmes robustes aux erreurs d'arrondis en abordant la comparaison de flottants, la résolution stable du second degré, la recherche de racines par dichotomie (bisection) et les méthodes classiques d'intégration numérique par quadrature.

10.2 Égalité de Flottants et Tolérances

Comme vu au Chapitre 3, deux calculs mathématiquement identiques peuvent produire des flottants légèrement différents en machine (ex. 0.1 + 0.2 != 0.3). Il est donc extrêmement risqué de baser la terminaison d'une boucle ou un test de validité sur une égalité stricte x == y.

10.2.1 Erreur Absolue vs Erreur Relative

Pour définir si deux flottants $x$ et $y$ sont proches, on peut employer :

Une tolérance absolue : $∣ x - y ∣ \leq ε_{a}$ . Ce critère est inadapté si les nombres comparés sont très petits devant $ε_{a}$ ou au contraire extrêmement grands.
Une tolérance relative : $∣ x - y ∣ \leq ε_{r} ∣ y ∣$ . Ce critère résout le problème d'échelle mais présente deux défauts : il n'est pas symétrique (on peut avoir $∣ x - y ∣ \leq ε_{r} ∣ y ∣$ mais $∣ y - x ∣ > ε_{r} ∣ x ∣$ ) et il échoue lorsque $x$ et $y$ sont proches de 0 ou de signes opposés.

En pratique, on combine les deux approches :

Définition 10.1Fonctions isclose

numpy.isclose(x, y) : $∣ x - y ∣ \leq ε_{a} + ε_{r} ∣ y ∣$ (par défaut $ε_{r} = 1 0^{- 5}$ , $ε_{a} = 1 0^{- 8}$ )
math.isclose(x, y) (Python 3.5+) : $∣ x - y ∣ \leq ε_{a} + ε_{r} max (∣ x ∣, ∣ y ∣)$ (par défaut $ε_{r} = 1 0^{- 9}$ , $ε_{a} = 0.0$ )

Exemple 10.2Utilisation de isclose


import math

print(0.1 + 0.2 == 0.3) # False
print(math.isclose(0.1 + 0.2, 0.3)) # True (tolérance par défaut)
print(math.isclose(0.1 + 0.2, 0.3, rel_tol=1e-10, abs_tol=1e-12))

10.3 Application : Résolution Stable du Second Degré

Le problème de l'égalité approchée se pose lors du calcul du discriminant $Δ = b^{2} - 4 a c$ d'une équation $a x^{2} + b x + c = 0$ .

10.3.1 Le Problème de la Cancellation

Soit l'équation $0.01 x^{2} + 0.2 x + 1 = 0$ (où le discriminant théorique est nul, $Δ = 0. 2^{2} - 4 \times 0.01 \times 1 = 0$ , admettant une racine double $x = - 10.0$ ).

L'évaluation naïve donne $Δ \approx 6.93 \cdot 1 0^{- 18} > 0$ , conduisant à deux racines simples erronées.
Pour l'équation $0.011025 x^{2} + 0.21 x + 1 = 0$ ( $Δ = 0$ théoriquement), la machine calcule $Δ \approx - 6.93 \cdot 1 0^{- 18} < 0$ , concluant faussement à l'absence de solution réelle.

10.3.2 Solution Numérique Stable

Pour corriger ce défaut, nous devons tester la proximité à zéro de $Δ$ relativement au terme $b^{2}$ . En posant une tolérance égale à l'epsilon machine ( $ε = 2^{- 52} \approx 2.22 \cdot 1 0^{- 16}$ ), la condition d'égalité à zéro devient $∣Δ∣ \leq ε b^{2}$ :

Exemple 10.3Résolution quadratique stable


import math

def solve_quadratique_stable(a, b, c, eps=2**(-52)):
    """
    Resout ax^2 + bx + c = 0 de maniere numeriquement stable.
    """
    if abs(a) < 1e-15: # Cas dégénéré (équation linéaire)
        if abs(b) < 1e-15:
            return None # Pas d'équation
        return [-c / b]
        
    delta = b * b - 4 * a * c
    seuil = eps * b * b # Tolérance relative à b²
    
    if delta < -seuil:
        return [] # Pas de racine réelle
    elif delta > seuil:
        sqrt_delta = math.sqrt(delta)
        
        # Formule stable (évite la cancellation pour la grande racine)
        if b >= 0:
            x1 = (-b - sqrt_delta) / (2 * a)
            x2 = (2 * c) / (-b - sqrt_delta)
        else:
            x1 = (-b + sqrt_delta) / (2 * a)
            x2 = (2 * c) / (-b + sqrt_delta)
        return sorted([x1, x2])
    else:
        # Racine double (discriminant proche de zéro)
        return [-b / (2 * a)]

# Test
print(solve_quadratique_stable(0.01, 0.2, 1)) # [-10.0]
print(solve_quadratique_stable(1, -3, 2)) # [1.0, 2.0]

10.4 Recherche de Racine par Dichotomie

La méthode de dichotomie (bisection) est un algorithme robuste permettant de trouver une solution approchée à l'équation $f (x) = 0$ pour une fonction continue $f$ sur un intervalle $[a, b]$ .

10.4.1 Principe et Algorithme

Le théorème des valeurs intermédiaires (TVI) assure que si $f (a)$ et $f (b)$ sont de signes opposés ( $f (a) f (b) \leq 0$ ), alors il existe au moins une racine $c \in [a, b]$ telle que $f (c) = 0$ .

L'algorithme de dichotomie construit deux suites de bornes $(u_{n})$ et $(v_{n})$ telles que :

Initialisation : $u_{0} = a$ et $v_{0} = b$ .
À chaque étape, on calcule le point milieu $m = \frac{u _{n} + v _{n}}{2}$ .
Si $f (u_{n}) f (m) \leq 0$ , la racine se trouve dans la moitié gauche : $u_{n + 1} = u_{n}$ , $v_{n + 1} = m$ .
Sinon, elle est dans la moitié droite : $u_{n + 1} = m$ , $v_{n + 1} = v_{n}$ .

10.4.2 Analyse de l'Erreur et Convergence

La largeur de l'intervalle est divisée par 2 à chaque étape :

$v_{n} - u_{n} = \frac{b - a}{2 ^{n}}$

L'invariant de boucle est $f (u_{n}) f (v_{n}) \leq 0$ , ce qui garantit la présence d'une racine dans $[u_{n}, v_{n}]$ .

Si l'on s'arrête lorsque $v_{n} - u_{n} \leq 2 ε$ , alors l'approximation finale $w_{n} = \frac{u _{n} + v _{n}}{2}$ vérifie $∣ w_{n} - c ∣ \leq ε$ .

Le nombre d'itérations $n$ nécessaire pour atteindre cette précision satisfait :

$\frac{b - a}{2 ^{n}} \leq 2 ε ⟺ n \geq lo g_{2} (\frac{b - a}{ε}) - 1$

10.4.3 Implémentation Python Sécurisée

Exemple 10.4Algorithme de dichotomie


def dichotomie(f, a, b, epsilon=1e-12, maxiter=100):
    """
    Recherche une racine de f par dichotomie sur [a, b].
    
    Paramètres:
    f : fonction continue
    a, b : bornes de l'intervalle (f(a)*f(b) <= 0)
    epsilon : précision souhaitée
    maxiter : nombre maximal d'itérations
    
    Retourne:
    Approximation de la racine
    """
    if f(a) * f(b) > 0:
        raise ValueError("f(a) et f(b) doivent être de signes opposés")
        
    u, v = a, b
    for n in range(maxiter):
        w = (u + v) / 2
        if abs(v - u) <= 2 * epsilon:
            return w
            
        if f(u) * f(w) <= 0:
            v = w
        else:
            u = w
            
    raise RuntimeError(f"Non convergence après {maxiter} itérations")

# Exemple : racine de x^2 - 2 (sqrt(2)) sur [1, 2]
f = lambda x: x**2 - 2
racine = dichotomie(f, 1, 2, epsilon=1e-12)
print(f"sqrt(2) = {racine:.12f}")

10.4.4 Méthode de Newton (Convergence Quadratique)

Lorsque la dérivée $f^{'}$ est disponible, la méthode de Newton offre une convergence quadratique.

Exemple 10.5Méthode de Newton


def newton(f, df, x0, epsilon=1e-12, maxiter=100):
    """
    Recherche une racine de f par la méthode de Newton.
    
    Paramètres:
    f : fonction
    df : dérivée de f
    x0 : approximation initiale
    epsilon : précision souhaitée
    maxiter : nombre maximal d'itérations
    """
    x = x0
    for n in range(maxiter):
        fx = f(x)
        if abs(fx) < epsilon:
            return x
        dfx = df(x)
        if abs(dfx) < 1e-15:
            raise RuntimeError("Dérivée nulle")
        x = x - fx / dfx
    raise RuntimeError(f"Non convergence après {maxiter} itérations")

# Exemple : sqrt(2) par Newton
f = lambda x: x**2 - 2
df = lambda x: 2*x
racine_newton = newton(f, df, 1.5)
print(f"sqrt(2) (Newton) = {racine_newton:.12f}")

10.4.5 Utilisation de SciPy

Exemple 10.6Recherche de racines avec SciPy


from scipy.optimize import bisect, newton

# Dichotomie SciPy
racine_bisect = bisect(lambda x: x**2 - 2, 1, 2, xtol=1e-12)

# Newton SciPy
racine_newton_sp = newton(lambda x: x**2 - 2, 1.5)

print(f"SciPy bisect : {racine_bisect:.12f}")
print(f"SciPy newton : {racine_newton_sp:.12f}")

10.5 Intégration Numérique (Méthodes de Quadrature)

Le calcul analytique d'une intégrale définie $I (f) = \int_{a}^{b} f (t) d t$ est parfois complexe ou impossible. Les méthodes de quadrature approchent cette intégrale par une somme pondérée finie de valeurs de la fonction :

$I_{u, v} (f) \approx i = 0 \sum p α_{i} f (x_{i})$

Pour calculer l'intégrale sur un grand intervalle $[a, b]$ , on le subdivise en $n$ sous-intervalles de taille $h = \frac{b - a}{n}$ délimités par les points $x_{k} = a + k h$ (ou $u_{k} = a + k h$ selon la notation).

10.5.1 Méthode du Rectangle (Gauche et Droite)

La méthode du rectangle à gauche sur un intervalle élémentaire $[u, v]$ approxime la fonction par sa valeur sur la borne inférieure :

$I_{u, v}^{rect-gauche} (f) = (v - u) f (u)$

La méthode du rectangle à droite utilise la borne supérieure :

$I_{u, v}^{rect-droite} (f) = (v - u) f (v)$

◆Théorème 10.7Erreur du Rectangle Élémentaire

Si $f$ est de classe $C^{1}$ sur $[u, v]$ , alors l'erreur vérifie :

$∣ E_{u, v} (f) ∣ \leq M_{1} \frac{( v - u ) ^{2}}{2}$

où $M_{1} = max_{[u, v]} ∣ f^{'} ∣$ .

Formules composites :

Rectangle à gauche :

$\int_{a}^{b} f (t) d t \approx h k = 0 \sum n - 1 f (x_{k})$

où $x_{k} = a + k h$ . L'erreur globale est en $O (1/ n)$ : $∣ E_{n} (f) ∣ \leq M_{1} \frac{( b - a ) ^{2}}{2 n}$ .

Rectangle à droite :

$\int_{a}^{b} f (t) d t \approx h k = 1 \sum n f (x_{k})$

Exemple 10.8Rectangle gauche et droit


def rectangle_gauche(f, a, b, n):
    h = (b - a) / n
    x = np.linspace(a, b - h, n)
    return h * np.sum(f(x))

def rectangle_droit(f, a, b, n):
    h = (b - a) / n
    x = np.linspace(a + h, b, n)
    return h * np.sum(f(x))

10.5.2 Méthode du Point Milieu

La méthode du point milieu approche l'intégrale en évaluant la fonction au centre de l'intervalle $w = \frac{u + v}{2}$ :

$I_{u, v}^{P M} (f) = (v - u) f (w)$

◆Théorème 10.9Erreur du Point Milieu

Si $f$ est de classe $C^{2}$ sur $[u, v]$ , alors l'erreur vérifie :

$∣ E_{u, v} (f) ∣ \leq M_{2} \frac{( v - u ) ^{3}}{24}$

où $M_{2} = max_{[u, v]} ∣ f^{''} ∣$ . La méthode est d'ordre 2.

Formule composite :

$\int_{a}^{b} f (t) d t \approx h k = 0 \sum n - 1 f (a + (k + \frac{1}{2}) h)$

L'erreur globale est en $O (1/ n^{2})$ : $∣ E_{n} (f) ∣ \leq M_{2} \frac{( b - a ) ^{3}}{24 n ^{2}}$ .

Exemple 10.10Méthode du point milieu


def point_milieu(f, a, b, n):
    h = (b - a) / n
    x = np.linspace(a + h/2, b - h/2, n)
    return h * np.sum(f(x))

10.5.3 Méthode des Trapèzes

La méthode des trapèzes approxime la fonction par la fonction affine joignant $(u, f (u))$ et $(v, f (v))$ :

$I_{u, v}^{t r a p} (f) = (v - u) \frac{f ( u ) + f ( v )}{2}$

◆Théorème 10.11Erreur du Trapèze Élémentaire

Si $f$ est de classe $C^{2}$ sur $[u, v]$ , alors l'erreur vérifie :

$∣ E_{u, v} (f) ∣ \leq M_{2} \frac{( v - u ) ^{3}}{12}$

La méthode est d'ordre 2.

Formule composite :

$\int_{a}^{b} f (t) d t \approx h (\frac{f ( a ) + f ( b )}{2} + k = 1 \sum n - 1 f (x_{k}))$

L'erreur globale est en $O (1/ n^{2})$ : $∣ E_{n} (f) ∣ \leq M_{2} \frac{( b - a ) ^{3}}{12 n ^{2}}$ .

Exemple 10.12Méthode des trapèzes


def trapezes(f, a, b, n):
    h = (b - a) / n
    x = np.linspace(a, b, n+1)
    return h * (np.sum(f(x)) - (f(a) + f(b))/2)

10.5.4 Méthode de Simpson

La méthode de Simpson approche la fonction par le polynôme de degré 2 passant par $(u, f (u))$ , $(w, f (w))$ et $(v, f (v))$ :

$I_{u, v}^{Simpson} (f) = \frac{v - u}{6} (f (u) + 4 f (w) + f (v))$

◆Théorème 10.13Erreur de Simpson

Bien que basée sur une interpolation quadratique (degré 2), la méthode de Simpson est d'ordre 4 (exacte pour les polynômes de degré $\leq 3$ ). Si $f$ est de classe $C^{4}$ sur $[u, v]$ , alors :

$∣ E_{u, v} (f) ∣ \leq M_{4} \frac{( v - u ) ^{5}}{2880}$

où $M_{4} = max_{[u, v]} ∣ f^{(4)} ∣$ .

Formule composite ( $n$ doit être pair) :

$\int_{a}^{b} f (t) d t \approx \frac{h}{3} f (a) + f (b) + 4 k = 0 \sum n /2 - 1 f (x_{2 k + 1}) + 2 k = 1 \sum n /2 - 1 f (x_{2 k})$

L'erreur globale est en $O (1/ n^{4})$ : $∣ E_{n} (f) ∣ \leq M_{4} \frac{( b - a ) ^{5}}{2880 n ^{4}}$ .

La suite du cours est sur Adloun

Le chapitre en entier, formules et figures comprises, avec le tuteur, les fiches, les QCM et les exercices.

Lire la suite du cours

En accès libre, sans inscription.

Exercices corrigés de ce chapitre

Lire la suite du cours →