Recherche de zéro d’une fonction

1Dichotomie et methode de Newton

Methode de dichotomie

def dichotomie(f, a, b, eps=1e-12, nmax=200):
    """Zero de f dans [a,b]. Precondition : f continue et f(a)*f(b) < 0.

    Renvoie (racine, nombre d'iterations).
    """
    fa, fb = f(a), f(b)
    if fa == 0: return a, 0
    if fb == 0: return b, 0
    if fa * fb > 0:
        raise ValueError("f(a) et f(b) doivent etre de signes opposes")
    n = 0
    while b - a > eps and n < nmax:
        m = (a + b) / 2
        fm = f(m)
        if fm == 0:
            return m, n
        if fa * fm < 0:
            b = m                      # le zero est dans [a, m]
        else:
            a, fa = m, fm              # il est dans [m, b]
        n += 1
    return (a + b) / 2, n

from math import cos
r, n = dichotomie(lambda x: cos(x) - x, 0, 1)
print(f"{r:.15f} en {n} iterations")     # 0.739085133215161 en 40 iterations

Methode de Newton

def newton(f, fp, x0, eps=1e-14, nmax=100):
    """Zero de f par la methode de Newton. fp est la derivee.

    Renvoie (racine, nombre d'iterations, convergence).
    """
    x = x0
    for n in range(1, nmax + 1):
        d = fp(x)
        if d == 0:
            return x, n, False           # tangente horizontale : echec
        x_nouveau = x - f(x) / d
        if abs(x_nouveau - x) < eps * max(1.0, abs(x_nouveau)):
            return x_nouveau, n, True
        x = x_nouveau
    return x, nmax, False

from math import cos, sin
r, n, ok = newton(lambda x: cos(x) - x, lambda x: -sin(x) - 1, 1.0)
print(f"{r:.15f} en {n} iterations")     # 0.739085133215161 en 5 iterations

from math import atan
# 1. divergence : le point de depart est trop loin
r, n, ok = newton(atan, lambda x: 1/(1+x*x), 2.0)
print(ok)          # False : la suite oscille en s'eloignant

# 2. tangente horizontale
newton(lambda x: x**3 - 2*x + 2, lambda x: 3*x*x - 2, 0.0)   # cycle 0 -> 1 -> 0

# 3. racine double : convergence seulement lineaire
newton(lambda x: (x-1)**2, lambda x: 2*(x-1), 2.0)

2Visualisation, secante et comparaison

Methode de la secante

Newton exige de connaitre \(f'\). Quand elle est inconnue ou couteuse, on la remplace par le taux d'accroissement entre les deux derniers points :

\[ x_{n+1} = x_n - f(x_n)\,\frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}. \]

def secante(f, x0, x1, eps=1e-14, nmax=100):
    """Newton sans derivee. Ordre de convergence : le nombre d'or ~ 1.618."""
    f0, f1 = f(x0), f(x1)
    for n in range(1, nmax + 1):
        if f1 == f0:
            return x1, n, False
        x2 = x1 - f1 * (x1 - x0) / (f1 - f0)
        if abs(x2 - x1) < eps * max(1.0, abs(x2)):
            return x2, n, True
        x0, f0, x1, f1 = x1, f1, x2, f(x2)
    return x1, nmax, False

Comparaison

3Methodes hybrides, criteres d'arret et applications

Methodes hybrides

def hybride(f, fp, a, b, eps=1e-14, nmax=100):
    """Newton securise par dichotomie. Convergence garantie ET rapide."""
    fa, fb = f(a), f(b)
    if fa * fb > 0:
        raise ValueError("pas de changement de signe")
    x = (a + b) / 2
    for n in range(1, nmax + 1):
        d = fp(x)
        pas_newton = (d != 0)
        if pas_newton:
            xn = x - f(x) / d
            pas_newton = (a < xn < b)          # le pas reste-t-il dans l'encadrement ?
        x = xn if pas_newton else (a + b) / 2   # sinon : dichotomie
        fx = f(x)
        if fx == 0 or (b - a) < eps:
            return x, n, "Newton" if pas_newton else "dichotomie"
        if fa * fx < 0:
            b = x
        else:
            a, fa = x, fx
    return x, nmax, "max atteint"

from math import atan
# atan fait diverger Newton depuis x0 = 2, mais l'hybride converge
print(hybride(atan, lambda x: 1/(1+x*x), -6, 6))

C'est le principe de la methode de Brent, utilisee par scipy.optimize.brentq : elle combine dichotomie, secante et interpolation quadratique inverse, avec convergence garantie et un ordre pratique proche de 1,8.

Choix du critere d'arret

while abs(f(x)) > eps:        # 1. residu : trompeur si |f'| est tres petit
while abs(x1 - x0) > eps:     # 2. increment absolu : mal adapte aux grandes valeurs
while abs(x1 - x0) > eps * abs(x1):   # 3. increment relatif : le plus robuste

def critere(x_nouveau, x, f_x, eps=1e-14, tol=1e-14):
    return (abs(x_nouveau - x) < eps * max(1.0, abs(x_nouveau))
            or abs(f_x) < tol)

Applications classiques

def racine_kieme(A, k, eps=1e-15):
    """Racine k-ieme de A > 0 par Newton (methode de Heron generalisee)."""
    x = A if A > 1 else 1.0
    for _ in range(100):
        x_nouveau = ((k-1)*x + A / x**(k-1)) / k
        if abs(x_nouveau - x) < eps * x:
            return x_nouveau
        x = x_nouveau
    return x

assert abs(racine_kieme(2, 2) - 2**0.5) < 1e-14
assert abs(racine_kieme(1000, 3) - 10) < 1e-12

def inverse(A, x0=None, n=30):
    """1/A par Newton, sans aucune division. Convergence si 0 < x0 < 2/A."""
    x = x0 if x0 is not None else 2.0**(-A.bit_length()) if isinstance(A, int) else 0.1
    for _ in range(n):
        x = x * (2 - A*x)
    return x

print(inverse(7, 0.1), 1/7)      # 0.14285714285714285 0.14285714285714285

def point_fixe(g, x0, eps=1e-14, nmax=1000):
    """Iteration x -> g(x). Converge si g est contractante au voisinage du point fixe."""
    x = x0
    for n in range(1, nmax+1):
        y = g(x)
        if abs(y - x) < eps * max(1.0, abs(y)):
            return y, n, True
        x = y
    return x, nmax, False

from math import cos, sqrt
print(point_fixe(cos, 1.0))                     # 0.739085..., ~85 iterations
print(point_fixe(lambda x: sqrt(2 + x), 0.0))   # 2.0, ~30 iterations