Méthodes numériques en Python : dichotomie, Newton, trapèzes

Quand une équation f(x) = 0 ou une intégrale ne se calcule pas à la main, on en cherche une valeur approchée par un algorithme : balayage, dichotomie ou méthode de Newton pour l'équation ; point milieu, trapèzes ou Monte-Carlo pour l'intégrale. Le programme du BTS cite ces méthodes et demande d'en programmer au moins une. Ce chapitre les code en Python, sur deux exemples qu'aucune formule ne résout, et compare leur précision.

Deux problèmes sans solution exacte

Une équation. Résoudre eˣ + x − 3 = 0. Aucune manipulation avec ln ne sort x de l'exponentielle et du terme x à la fois. Mais on sait qu'il y a une solution unique : f(x) = eˣ + x − 3 a pour dérivée eˣ + 1 > 0, elle est strictement croissante, et f(0) = −2 < 0 < f(1) = e − 2 ≈ 0,718. Il y a donc une seule solution α, entre 0 et 1 (chapitre 3).

Une intégrale. Calculer ∫₀¹ e^(−x²) dx. La fonction e^(−x²) est continue et positive, l'aire sous sa courbe existe, mais elle n'a aucune primitive qui s'écrive avec les fonctions usuelles. À un changement d'échelle près, c'est la courbe en cloche de la loi normale : c'est pour cela que la calculatrice, et non une formule, calcule les probabilités de la loi normale.

Le balayage

On avance de x = a à x = b par pas réguliers, et on s'arrête au premier intervalle où f change de signe, c'est-à-dire où f(x) × f(x + pas) ≤ 0.

from math import exp


def f(x):
    return exp(x) + x - 3


def balayage(f, a, b, pas):
    """Premier intervalle [x ; x + pas] de [a ; b] sur lequel f change de signe."""
    n = round((b - a) / pas)          # nombre de tranches
    for i in range(n):
        x = round(a + i * pas, 10)    # arrondi : 0.1 n'est pas exact en binaire
        if f(x) * f(x + pas) <= 0:
            return x, round(x + pas, 10)
    return None


print(balayage(f, 0, 1, 0.1))
print(balayage(f, 0.7, 0.8, 0.01))

(0.7, 0.8)
(0.79, 0.8)

Premier passage au pas de 0,1 : α est entre 0,7 et 0,8. Deuxième passage au pas de 0,01 : α est entre 0,79 et 0,80. Chaque chiffre de plus coûte un nouveau balayage d'une dizaine de calculs.

La dichotomie

La dichotomie coupe l'intervalle en deux à chaque étape, comme la recherche dichotomique dans un tableau trié (chapitre 6) : on calcule f au milieu m, et on garde la moitié où f change de signe.

from math import exp


def f(x):
    return exp(x) + x - 3


def dichotomie(f, a, b, precision):
    """Encadrement [a ; b] d'amplitude au plus precision d'une solution de f(x) = 0.
    Précondition : f est continue sur [a ; b] et f(a), f(b) sont de signes contraires."""
    etapes = 0
    while b - a > precision:
        m = (a + b) / 2
        if f(a) * f(m) <= 0:   # le changement de signe est dans [a ; m]
            b = m
        else:                  # sinon, il est dans [m ; b]
            a = m
        etapes = etapes + 1
    return a, b, etapes


print(dichotomie(f, 0, 1, 0.001))
print(dichotomie(f, 0, 1, 1e-6))

(0.7919921875, 0.79296875, 10)
(0.7920598983764648, 0.7920608520507812, 20)

Avec 10 étapes, α est connu à 10⁻³ près ; avec 20 étapes, à 10⁻⁶ près : α ≈ 0,792 06.

La dichotomie est lente mais sûre : tant que f est continue et change de signe sur l'intervalle de départ, elle converge, et l'on sait d'avance combien d'étapes il faut.

La méthode de Newton

La méthode de Newton remplace la courbe par sa tangente. En partant d'une valeur x₀, on trace la tangente au point d'abscisse x₀ ; elle coupe l'axe des abscisses en x₁, plus proche de la solution ; on recommence. L'équation de la tangente (chapitre 2) donne la formule :

xₙ₊₁ = xₙ − f(xₙ) / f′(xₙ)

from math import exp


def f(x):
    return exp(x) + x - 3


def derivee_f(x):
    return exp(x) + 1


def newton(f, df, x0, n):
    """Calcule n itérations de la méthode de Newton à partir de x0."""
    x = x0
    for i in range(1, n + 1):
        x = x - f(x) / df(x)
        print(i, x)
    return x


newton(f, derivee_f, 1, 4)

1 0.8068242641099854
2 0.7921349596563866
3 0.7920599703659901
4 0.792059968430677

En trois itérations, Newton donne α à 10⁻⁸ près ; la dichotomie demande 27 étapes pour la même précision (log₂(10⁸) ≈ 26,6). En contrepartie, il faut connaître f′, et un mauvais point de départ (près d'un point où f′ s'annule) peut faire diverger la méthode.

Approcher une intégrale : rectangles et point milieu

Pour approcher ∫ₐᵇ f(x) dx, on découpe [a ; b] en n tranches de largeur h = (b − a)/n et on remplace l'aire sous la courbe, sur chaque tranche, par celle d'un rectangle. Dans la méthode du point milieu, la hauteur du rectangle est la valeur de f au milieu de la tranche :

∫ₐᵇ f(x) dx ≈ h × [f(a + h/2) + f(a + 3h/2) + … + f(b − h/2)]

La méthode des trapèzes

Dans la méthode des trapèzes, on relie les points de la courbe aux bords de chaque tranche par un segment : chaque tranche devient un trapèze, d'aire h × (f(gauche) + f(droite))/2. En additionnant, les valeurs intérieures comptent deux fois :

∫ₐᵇ f(x) dx ≈ h × [f(a)/2 + f(a + h) + f(a + 2h) + … + f(b − h) + f(b)/2]

La méthode de Monte-Carlo

La méthode de Monte-Carlo tire des points au hasard dans un rectangle qui contient le domaine, et compte la proportion qui tombe sous la courbe. Cette proportion, multipliée par l'aire du rectangle, approche l'intégrale. Elle est peu précise, mais elle se généralise facilement aux domaines compliqués et aux dimensions supérieures.

Programmer les trois méthodes en Python

Chaque méthode tient en quelques lignes. Le programme ci-dessous les compare sur ∫₀¹ e^(−x²) dx.

from math import exp
import random


def f(x):
    return exp(-x ** 2)


def point_milieu(f, a, b, n):
    """Somme des aires de n rectangles dont la hauteur est prise au milieu de chaque tranche."""
    h = (b - a) / n
    s = 0
    for i in range(n):
        s = s + f(a + (i + 0.5) * h)
    return s * h


def trapezes(f, a, b, n):
    """Somme des aires de n trapèzes."""
    h = (b - a) / n
    s = (f(a) + f(b)) / 2
    for i in range(1, n):
        s = s + f(a + i * h)
    return s * h


def monte_carlo(f, a, b, ymax, n):
    """Proportion de n points tirés au hasard dans le rectangle [a ; b] x [0 ; ymax]
    qui tombent sous la courbe, multipliée par l'aire du rectangle.
    Précondition : 0 <= f(x) <= ymax sur [a ; b]."""
    dessous = 0
    for _ in range(n):
        x = random.uniform(a, b)
        y = random.uniform(0, ymax)
        if y <= f(x):
            dessous = dessous + 1
    return dessous / n * (b - a) * ymax


for n in (4, 10, 100):
    print(n, point_milieu(f, 0, 1, n), trapezes(f, 0, 1, n))

random.seed(2026)
print(monte_carlo(f, 0, 1, 1, 100_000))

4 0.7487471318910093 0.7429840978003812
10 0.7471308777479975 0.7462107961317495
100 0.7468271984923199 0.7468180014679697
0.75041

La valeur exacte, au millionième, est 0,746 824. Comparons :

• Avec 100 tranches, le point milieu et les trapèzes ont cinq décimales justes. Multiplier n par 10 divise l'erreur par environ 100.
• Ici, le point milieu est un peu plus précis que les trapèzes, et les deux erreurs sont de signes contraires : la valeur exacte est encadrée.
• Monte-Carlo, avec 100 000 points, se trompe encore d'environ 0,004, et son résultat change à chaque exécution (la ligne random.seed(2026) fixe le tirage pour que tu obtiennes le même nombre que ce cours).

Exercices

def rectangles_gauche(f, a, b, n):
    h = (b - a) / n
    s = 0
    for i in range(n):
        s = s + f(a + i * h)   # bord gauche de la tranche i
    return s * h


print(rectangles_gauche(lambda x: x ** 2, 0, 2, 100))