Recherche de zéro d’une fonction

1. Dichotomie (\(\ell = \sqrt2 \approx 1{,}414214\)) :

2. Newton : \(x_{n+1} = x_n - \frac{x_n^2-2}{2x_n} = \frac12(x_n + 2/x_n)\).

3. Apres 4 iterations : dichotomie \(2{,}3\cdot10^{-2}\), Newton \(1{,}6\cdot10^{-12}\). Dix ordres de grandeur d'ecart. Le nombre de decimales exactes de Newton suit 0, 1, 2, 5, 11 : il double a chaque etape.

4.

• Dichotomie : \(|x_n - \ell| \le 2^{-(n+1)}\), d'ou \(2^{-(n+1)} \le 10^{-10}\), soit \(n \ge 10\log_2 10 - 1 \approx 32\) iterations.
• Newton : avec \(e_{n+1} \approx 0{,}354\,e_n^2\), on part de \(e_0 = 0{,}59\) et l'on obtient \(e_1 \approx 8{,}6\cdot10^{-2}\), \(e_2 \approx 2{,}5\cdot10^{-3}\), \(e_3 \approx 2{,}1\cdot10^{-6}\), \(e_4 \approx 1{,}6\cdot10^{-12}\) : 4 iterations suffisent.

Rapport : 8 fois moins d'iterations, chacune coutant 2 evaluations au lieu d'une, soit un gain net de facteur 4.

from math import cos, sin, exp, log

def dichotomie(f, a, b, eps=1e-14, nmax=200):
    fa = f(a); ev = 2
    if fa * f(b) > 0:
        raise ValueError("pas de changement de signe")
    n = 0
    while b - a > eps * max(1.0, abs(a)) and n < nmax:
        m = (a+b)/2; fm = f(m); ev += 1; n += 1
        if fm == 0: return m, n, ev
        if fa*fm < 0: b = m
        else: a, fa = m, fm
    return (a+b)/2, n, ev

def newton(f, fp, x0, eps=1e-14, nmax=100):
    x, ev = x0, 0
    for n in range(1, nmax+1):
        fx, d = f(x), fp(x); ev += 2
        if d == 0: return x, n, ev
        y = x - fx/d
        if abs(y-x) < eps*max(1.0, abs(y)): return y, n, ev
        x = y
    return x, nmax, ev

def secante(f, x0, x1, eps=1e-14, nmax=100):
    f0, f1 = f(x0), f(x1); ev = 2
    for n in range(1, nmax+1):
        if f1 == f0: return x1, n, ev
        x2 = x1 - f1*(x1-x0)/(f1-f0)
        if abs(x2-x1) < eps*max(1.0, abs(x2)): return x2, n, ev
        x0, f0 = x1, f1
        x1, f1 = x2, f(x2); ev += 1
    return x1, nmax, ev

problemes = [
    ("cos(x) - x",  lambda x: cos(x)-x,   lambda x: -sin(x)-1, 0, 1, 0.5, 1.0),
    ("x^3 - 2",     lambda x: x**3-2,     lambda x: 3*x*x,     1, 2, 1.0, 2.0),
    ("exp(x) - 3",  lambda x: exp(x)-3,   exp,                 0, 2, 1.0, 2.0),
    ("x ln x - 1",  lambda x: x*log(x)-1, lambda x: log(x)+1,  1, 3, 1.5, 2.0),
]

print(f"{'probleme':<14}{'methode':<12}{'racine':>20}{'iter':>6}{'eval f':>8}")
for nom, f, fp, a, b, x0, x1 in problemes:
    for meth, res in [("dichotomie", dichotomie(f, a, b)),
                      ("Newton",     newton(f, fp, x1)),
                      ("secante",    secante(f, x0, x1))]:
        r, n, ev = res
        print(f"{nom:<14}{meth:<12}{r:>20.15f}{n:>6}{ev:>8}")
    print()

probleme      methode                   racine  iter  eval f
cos(x) - x    dichotomie   0.739085133215161    47      49
cos(x) - x    Newton       0.739085133215161     5      10
cos(x) - x    secante      0.739085133215161     6       8

x^3 - 2       dichotomie   1.259921049894873    47      49
x^3 - 2       Newton       1.259921049894873     6      12
x^3 - 2       secante      1.259921049894873     8      10

exp(x) - 3    dichotomie   1.098612288668110    48      50
exp(x) - 3    Newton       1.098612288668110     6      12
exp(x) - 3    secante      1.098612288668110     8      10

x ln x - 1    dichotomie   1.763222834351896    48      50
x ln x - 1    Newton       1.763222834351896     5      10
x ln x - 1    secante      1.763222834351896     7       9

4. Conclusion.

• La dichotomie demande systematiquement une cinquantaine d'iterations, conformement a \(\log_2(1/10^{-14}) \approx 47\). Elle est previsible mais lente.
• Newton converge en 5 a 6 iterations, mais chaque iteration coute deux evaluations (\(f\) et \(f'\)).
• La secante demande un peu plus d'iterations mais moins d'evaluations totales : c'est la plus economique ici, et elle n'exige pas de connaitre \(f'\).

En pratique : secante ou Brent quand \(f'\) est inconnue ou couteuse ; Newton quand \(f'\) est gratuite (fonction polynomiale, differentiation automatique) ; dichotomie quand la robustesse prime absolument, ou pour fournir un bon point de depart aux autres.

from math import atan, tan

def newton_trace(f, fp, x0, n=10):
    x, hist = x0, [x0]
    for _ in range(n):
        d = fp(x)
        if d == 0: break
        x = x - f(x)/d
        hist.append(x)
        if abs(x) > 1e100: break
    return hist

# 1. arctan
for x0 in (1.0, 1.3, 1.39, 1.4, 1.5, 2.0):
    h = newton_trace(atan, lambda x: 1/(1+x*x), x0, 8)
    etat = "converge" if abs(h[-1]) < 1e-8 else "DIVERGE"
    print(f"x0 = {x0:>5} : {etat:<10} |x8| = {abs(h[-1]):.3e}")

x0 =   1.0 : converge   |x8| = 0.000e+00
x0 =   1.3 : converge   |x8| = 0.000e+00
x0 =  1.39 : converge   |x8| = 5.551e-17
x0 =   1.4 : DIVERGE    |x8| = 1.093e+04
x0 =   1.5 : DIVERGE    |x8| = 1.696e+08
x0 =   2.0 : DIVERGE    |x8| = 5.766e+62

1. L'iteration de Newton pour \(\arctan\) est \(x_{n+1} = x_n - (1+x_n^2)\arctan x_n\). Le seuil est la solution de \(2x = (1+x^2)\arctan x\), c'est-a-dire le point ou \(|x_1| = |x_0|\) : il vaut \(x^* \approx 1{,}3917\). Au-dela, \(|x_{n+1}| > |x_n|\) : la suite oscille en s'eloignant, et diverge tres vite (\(|x_n|\) croit comme une tour d'exponentielles).

from scipy.optimize import brentq   # ou une dichotomie a la main
seuil = brentq(lambda x: (1+x*x)*atan(x) - 2*x, 1.0, 2.0)
print(f"seuil = {seuil:.6f}")       # 1.391746

# 2. cycle
h = newton_trace(lambda x: x**3-2*x+2, lambda x: 3*x*x-2, 0.0, 8)
print([f"{v:.6f}" for v in h])      # ['0.000000', '1.000000', '0.000000', '1.000000', ...]

2. La suite entre dans un cycle de periode 2 : \(0 \to 1 \to 0 \to 1 \to \cdots\). En effet \(f(0) = 2\), \(f'(0) = -2\), donc \(x_1 = 1\) ; puis \(f(1) = 1\), \(f'(1) = 1\), donc \(x_2 = 0\). La methode ne converge jamais, sans jamais diverger non plus : aucun critere d'arret sur \(|x_{n+1}-x_n|\) ne se declenche.

# 3. racine double
h = newton_trace(lambda x: (x-1)**2, lambda x: 2*(x-1), 2.0, 12)
e = [abs(v-1) for v in h]
print("erreurs :", [f"{v:.3e}" for v in e[:6]])
print("rapports :", [f"{e[i+1]/e[i]:.4f}" for i in range(5)])

erreurs : ['1.000e+00', '5.000e-01', '2.500e-01', '1.250e-01', '6.250e-02', '3.125e-02']
rapports : ['0.5000', '0.5000', '0.5000', '0.5000', '0.5000']

3. Le rapport est constant egal a \(1/2\) : la convergence est lineaire, pas mieux que la dichotomie. Pour une racine de multiplicite \(m\), on montre que le rapport vaut \(1 - 1/m\).

4. Remedes.

def newton_modifie(f, fp, m, x0, n=20):
    x = x0
    for _ in range(n):
        x = x - m * f(x) / fp(x)
    return x

print(newton_modifie(lambda x: (x-1)**2, lambda x: 2*(x-1), 2, 2.0, 6))
# 1.0 : la convergence quadratique est retrouvee

from math import sqrt, log

exact = sqrt(2)

# 1-2. Newton
x, E = 2.0, []
for _ in range(6):
    x = 0.5*(x + 2/x)
    E.append(abs(x - exact))
print("erreurs :", [f"{v:.3e}" for v in E])
print("e(n+1)/e(n)^2 :", [f"{E[i+1]/E[i]**2:.6f}" for i in range(3)])
print(f"theorie : |f''|/(2|f'|) en sqrt(2) = {2/(2*2*sqrt(2)):.6f}")

erreurs : ['8.579e-02', '2.453e-03', '2.124e-06', '1.595e-12', '2.220e-16', '2.220e-16']
e(n+1)/e(n)^2 : ['0.333333', '0.352941', '0.353553']
theorie : |f''|/(2|f'|) en sqrt(2) = 0.353553

2. Pour \(f(x) = x^2-2\), on a \(f'(\ell) = 2\sqrt2\) et \(f'' = 2\), d'ou \[ K = \frac{|f''(\ell)|}{2|f'(\ell)|} = \frac{2}{4\sqrt2} = \frac{1}{2\sqrt2} = \frac{\sqrt2}{4} \approx 0{,}353553. \] Le rapport mesure converge exactement vers cette valeur.

# 3. secante
def secante_erreurs(f, x0, x1, n=8):
    f0, f1, E = f(x0), f(x1), []
    for _ in range(n):
        if f1 == f0: break
        x2 = x1 - f1*(x1-x0)/(f1-f0)
        E.append(abs(x2 - exact))
        x0, f0, x1, f1 = x1, f1, x2, f(x2)
    return E

E = secante_erreurs(lambda x: x*x-2, 1.0, 2.0)
print("erreurs :", [f"{v:.3e}" for v in E[:6]])
print("ordre estime :")
for i in range(1, 4):
    if E[i] > 1e-15 and E[i-1] > 1e-15:
        print(f"  {log(E[i+1]/E[i]) / log(E[i]/E[i-1]):.4f}")
print(f"nombre d'or : {(1+sqrt(5))/2:.4f}")

ordre estime :
  1.6154
  1.6180
  1.6180
nombre d'or : 1.6180

3. L'ordre converge vers \(\varphi = \frac{1+\sqrt5}{2}\). L'analyse : on montre \(e_{n+1} \approx \frac{|f''|}{2|f'|} e_n e_{n-1}\) ; en posant \(e_n \sim A e_{n-1}^\alpha\), il vient \(A e_{n-1}^{\alpha^2} \sim K A e_{n-1}^{\alpha+1}\), d'ou \(\alpha^2 = \alpha + 1\).

# 4. dichotomie : la LARGEUR est exactement divisee par 2
a, b, L = 1.0, 2.0, []
for _ in range(8):
    m = (a+b)/2
    if m*m > 2: b = m
    else: a = m
    L.append(b - a)
print("largeurs :", [f"{v:.6f}" for v in L])
print("rapports :", [f"{L[i+1]/L[i]:.6f}" for i in range(5)])   # tous egaux a 0.5

4. Le rapport des largeurs est exactement \(1/2\). L'erreur \(|m_n - \ell|\), elle, oscille (elle peut diminuer plus vite par chance a une etape donnee) mais reste majoree par la demi-largeur. C'est pourquoi on parle de convergence lineaire de rapport \(1/2\) en encadrement.

from math import sqrt, exp, log

def racine_carree(A, eps=1e-15):
    """sqrt(A) par la methode de Heron. A >= 0."""
    if A < 0: raise ValueError("A doit etre positif")
    if A == 0: return 0.0
    x = A if A >= 1 else 1.0            # initialisation raisonnable
    for _ in range(100):
        y = 0.5*(x + A/x)
        if abs(y - x) < eps*y: return y
        x = y
    return x

def racine_kieme(A, k, eps=1e-15):
    """A^(1/k) pour A > 0, k >= 1."""
    x = max(1.0, A**(1.0/k) if A < 1e300 else 1.0)
    for _ in range(200):
        y = ((k-1)*x + A/x**(k-1)) / k
        if abs(y - x) < eps*abs(y): return y
        x = y
    return x

def inverse(A, eps=1e-16, nmax=100):
    """1/A SANS division. Precondition : A != 0.

    Initialisation par mise a l'echelle : on ramene |A| dans [0.5, 1].
    """
    signe = 1.0 if A > 0 else -1.0
    A = abs(A)
    echelle = 1.0
    while A > 1.0:
        A *= 0.5; echelle *= 0.5
    while A < 0.5:
        A *= 2.0; echelle *= 2.0
    x = 2.914 - 2*A                     # bonne approximation initiale sur [0.5, 1]
    for _ in range(nmax):
        y = x * (2 - A*x)               # aucune division
        if abs(y - x) < eps: break
        x = y
    return signe * x * echelle

def log_newton(A, eps=1e-15):
    """ln(A) en resolvant exp(x) - A = 0. Newton : x -> x - 1 + A/exp(x).

    L'initialisation est CRUCIALE : partir de x0 = 0 pour A = 1e10 donnerait
    x1 = 1e10 - 1, puis exp(x1) = OverflowError. On utilise la decomposition
    A = m * 2^k (frexp) : ln(A) = ln(m) + k ln(2), avec m dans [0.5, 1[.
    """
    from math import frexp, log2
    if A <= 0:
        raise ValueError("A doit etre strictement positif")
    m, k = frexp(A)                     # A = m * 2^k, 0.5 <= m < 1
    x = (m - 1) + k * 0.6931471805599453    # ln(m) ~ m-1, plus k*ln(2)
    for _ in range(200):
        y = x - 1 + A / exp(x)          # = x - (exp(x) - A)/exp(x)
        if abs(y - x) < eps*max(1.0, abs(y)):
            return y
        x = y
    return x

for A in (2, 10, 1e-6, 1e10, 0.5):
    assert abs(racine_carree(A) - sqrt(A)) < 1e-12 * sqrt(A)
    assert abs(inverse(A) - 1/A) < 1e-12 / A
    assert abs(log_newton(A) - log(A)) < 1e-12 * max(1, abs(log(A)))
assert abs(racine_kieme(1000, 3) - 10) < 1e-12
assert abs(racine_kieme(2, 10) - 2**0.1) < 1e-14
print("toutes les fonctions validees")

from math import atan, cos, sin, exp

def hybride(f, fp, a, b, eps=1e-14, nmax=100):
    """Newton securise par encadrement. Renvoie (racine, nb_newton, nb_dicho)."""
    fa, fb = f(a), f(b)
    if fa*fb > 0:
        raise ValueError("pas de changement de signe sur [a, b]")
    x = (a + b)/2
    nb_n = nb_d = 0
    for _ in range(nmax):
        d = fp(x)
        candidat = None
        if d != 0:
            xn = x - f(x)/d
            if a < xn < b:
                candidat = xn
        if candidat is not None:
            x, nb_n = candidat, nb_n + 1
        else:
            x, nb_d = (a + b)/2, nb_d + 1
        fx = f(x)
        if fx == 0 or (b - a) < eps*max(1.0, abs(x)):
            return x, nb_n, nb_d
        if fa*fx < 0:
            b = x
        else:
            a, fa = x, fx
    return x, nb_n, nb_d

def newton(f, fp, x0, nmax=60):
    x = x0
    for _ in range(nmax):
        d = fp(x)
        if d == 0 or abs(x) > 1e100: return x, False
        x = x - f(x)/d
    return x, abs(f(x)) < 1e-10

datan = lambda x: 1/(1+x*x)
print("Newton seul depuis 6 :", newton(atan, datan, 6.0))          # diverge
r, nn, nd = hybride(atan, datan, -6, 6)
print(f"hybride : racine {r:.15f}, {nn} pas de Newton, {nd} pas de dichotomie")

for nom, f, fp, a, b in [("cos(x)-x", lambda x: cos(x)-x, lambda x: -sin(x)-1, -5, 5),
                         ("exp(x)-3", lambda x: exp(x)-3, exp, -10, 10),
                         ("x^3-2",    lambda x: x**3-2,   lambda x: 3*x*x, 0.1, 100)]:
    r, nn, nd = hybride(f, fp, a, b)
    print(f"{nom:<12} racine {r:.14f}, Newton {nn}, dichotomie {nd}")

Newton seul depuis 6 : (nan, False)
hybride : racine 0.000000000000000, 4 pas de Newton, 6 pas de dichotomie
cos(x)-x     racine 0.73908513321516, Newton 6, dichotomie 3
exp(x)-3     racine 1.09861228866811, Newton 7, dichotomie 5
x^3-2        racine 1.25992104989487, Newton 9, dichotomie 4

3. Le comportement typique : quelques pas de dichotomie au debut, tant que l'on est loin de la racine et que Newton propose des pas aberrants, puis Newton prend le relais et converge quadratiquement. On obtient a la fois la garantie de la dichotomie et la vitesse de Newton.

4. Comparaison avec Brent.

from scipy.optimize import brentq
import time

for nom, f, fp, a, b in [("cos(x)-x", lambda x: cos(x)-x, lambda x: -sin(x)-1, -5, 5)]:
    r1 = hybride(f, fp, a, b)[0]
    r2 = brentq(f, a, b, xtol=1e-15)
    print(f"{nom} : hybride {r1:.15f}, brentq {r2:.15f}, ecart {abs(r1-r2):.2e}")

La methode de Brent (1973) va plus loin : elle combine dichotomie, secante et interpolation quadratique inverse, choisit a chaque etape la plus prometteuse, et garantit un nombre d'iterations au plus egal au carre de celui de la dichotomie tout en atteignant en pratique un ordre proche de 1,8 — sans jamais avoir besoin de \(f'\). C'est l'algorithme par defaut de la plupart des bibliotheques scientifiques.

from math import sqrt

phi = (1 + sqrt(5)) / 2       # 1.618033988749895

def point_fixe(g, x0, nmax=60):
    x, hist = x0, [x0]
    for _ in range(nmax):
        try:
            x = g(x)
        except (ValueError, ZeroDivisionError, OverflowError):
            return hist, False
        hist.append(x)
        if abs(x) > 1e10: return hist, False
        if abs(x - phi) < 1e-14: return hist, True
    return hist, abs(hist[-1] - phi) < 1e-8

for nom, g, gp in [("g1(x) = x^2 - 1",   lambda x: x*x - 1,   lambda x: 2*x),
                   ("g2(x) = 1 + 1/x",   lambda x: 1 + 1/x,   lambda x: -1/(x*x)),
                   ("g3(x) = sqrt(1+x)", lambda x: sqrt(1+x), lambda x: 1/(2*sqrt(1+x)))]:
    h, ok = point_fixe(g, 1.5)
    print(f"{nom:<20} |g'(phi)| = {abs(gp(phi)):.4f}  "
          f"{'converge en ' + str(len(h)-1) + ' iterations' if ok else 'DIVERGE'}")

g1(x) = x^2 - 1      |g'(phi)| = 3.2361  DIVERGE
g2(x) = 1 + 1/x      |g'(phi)| = 0.3820  converge en 34 iterations
g3(x) = sqrt(1+x)    |g'(phi)| = 0.3090  converge en 30 iterations

2. Critere de convergence. Les trois fonctions ont bien \(\varphi\) pour point fixe (on le verifie en resolvant \(g(x) = x\)). Le theoreme du point fixe dit que l'iteration converge localement si et seulement si \(|g'(\ell)| < 1\) :

• \(g_1'(\varphi) = 2\varphi \approx 3{,}24 > 1\) : repulsif, l'iteration s'eloigne ;
• \(g_2'(\varphi) = -1/\varphi^2 \approx -0{,}382\) : attractif ;
• \(g_3'(\varphi) = \frac{1}{2\sqrt{1+\varphi}} = \frac{1}{2\varphi} \approx 0{,}309\) : attractif, et le plus rapide des trois puisque \(|g_3'| < |g_2'|\).

Le meme probleme, reecrit differemment, donne donc un algorithme qui converge ou qui diverge : la reformulation n'est pas neutre.

# 3. verification : e(n+1)/e(n) -> |g'(phi)|
for nom, g, gp in [("g2", lambda x: 1 + 1/x, lambda x: -1/(x*x)),
                   ("g3", lambda x: sqrt(1+x), lambda x: 1/(2*sqrt(1+x)))]:
    h, _ = point_fixe(g, 1.5)
    E = [abs(v - phi) for v in h if abs(v - phi) > 1e-15]
    r = [E[i+1]/E[i] for i in range(len(E)-1)]
    print(f"{nom} : rapports {[f'{v:.4f}' for v in r[-4:]]}, "
          f"|g'(phi)| = {abs(gp(phi)):.4f}")

g2 : rapports ['0.3820', '0.3820', '0.3820', '0.3820'], |g'(phi)| = 0.3820
g3 : rapports ['0.3090', '0.3090', '0.3090', '0.3090'], |g'(phi)| = 0.3090

Le rapport des erreurs converge exactement vers \(|g'(\ell)|\), comme le predit le developpement \(e_{n+1} = |g(x_n)-g(\ell)| \approx |g'(\ell)|\,e_n\).

4. Newton comme point fixe. Posons \(g(x) = x - \dfrac{f(x)}{f'(x)}\). Alors

\[ g'(x) = 1 - \frac{f'(x)^2 - f(x)f''(x)}{f'(x)^2} = \frac{f(x)\,f''(x)}{f'(x)^2}. \]

En \(\ell\), zero simple de \(f\), on a \(f(\ell) = 0\) et \(f'(\ell) \ne 0\), donc \(g'(\ell) = 0\).

Le facteur de contraction est nul : le terme lineaire du developpement de l'erreur disparait, et c'est le terme d'ordre 2 qui domine. On retrouve \[ e_{n+1} \approx \frac{|g''(\ell)|}{2}e_n^2 = \frac{|f''(\ell)|}{2|f'(\ell)|}e_n^2, \] c'est-a-dire exactement la convergence quadratique demontree en cours. Newton est le point fixe optimal : c'est le choix de \(g\) qui annule \(g'(\ell)\).

Si \(\ell\) est une racine double, alors \(f'(\ell) = 0\) aussi et l'on montre \(g'(\ell) = 1/2\) : on retombe sur une convergence lineaire de rapport \(1/2\), conformement a l'exercice 3.

1. Les racines de \(z^3 = 1\) sont \(1,\ j = e^{2\mathrm i\pi/3} = -\frac12 + \frac{\sqrt3}{2}\mathrm i,\ \bar j\).

2. L'iteration de Newton s'ecrit \[ z_{n+1} = z_n - \frac{z_n^3 - 1}{3z_n^2} = \frac{2z_n}{3} + \frac{1}{3z_n^2}. \] Elle est valable telle quelle dans \(\mathbb C\), la demonstration de la convergence quadratique etant identique (holomorphie au lieu de derivabilite reelle).

import numpy as np
import matplotlib.pyplot as plt

def newton_fractale(n=600, xmin=-2, xmax=2, ymin=-2, ymax=2, nmax=40, tol=1e-8):
    """Bassins d'attraction de Newton pour z^3 - 1. Entierement vectorise."""
    x = np.linspace(xmin, xmax, n)
    y = np.linspace(ymin, ymax, n)
    Z = x[None, :] + 1j * y[:, None]

    racines = np.array([1, np.exp(2j*np.pi/3), np.exp(-2j*np.pi/3)])
    iterations = np.zeros(Z.shape, dtype=int)
    actif = np.ones(Z.shape, dtype=bool)

    for k in range(nmax):
        Zk = Z[actif]
        Zk = Zk - (Zk**3 - 1) / (3*Zk**2 + 1e-300)   # on evite la division par 0
        Z[actif] = Zk
        proche = np.min(np.abs(Z[:, :, None] - racines[None, None, :]), axis=2) < tol
        nouveaux = actif & proche
        iterations[nouveaux] = k
        actif &= ~proche
        if not actif.any():
            break

    bassin = np.argmin(np.abs(Z[:, :, None] - racines[None, None, :]), axis=2)
    return bassin, iterations

bassin, it = newton_fractale(600)

fig, ax = plt.subplots(1, 2, figsize=(13, 6))
ax[0].imshow(bassin, extent=[-2, 2, -2, 2], cmap="Set1", origin="lower")
ax[0].set_title("Bassins d'attraction (une couleur par racine)")
# nuance par la vitesse de convergence
ax[1].imshow(bassin + 0.85*np.tanh(it/12), extent=[-2, 2, -2, 2],
             cmap="twilight", origin="lower")
ax[1].set_title("Bassins nuances par le nombre d'iterations")
for a in ax:
    a.set_xlabel("Re(z)"); a.set_ylabel("Im(z)")
fig.tight_layout(); plt.show()

# verification : les racines elles-memes convergent immediatement
for r in (1, np.exp(2j*np.pi/3), np.exp(-2j*np.pi/3)):
    assert abs(r**3 - 1) < 1e-14

4. Structure obtenue. Le plan se partage en trois bassins d'attraction, un par racine. Chacun est un ouvert dense au voisinage de sa racine, mais la frontiere entre les bassins est une fractale : l'ensemble de Julia de l'application de Newton.

La propriete remarquable est la suivante :

Tout point de la frontiere est adherent aux TROIS bassins simultanement. Il n'existe aucun arc de courbe separant seulement deux bassins.

C'est un resultat profond de la theorie de Fatou et Julia : la frontiere commune de trois ouverts disjoints ne peut etre une courbe reguliere (analogue au « lac de Wada »). Consequence pratique : au voisinage de la frontiere, une perturbation arbitrairement petite du point de depart change la racine atteinte. Le probleme « vers quelle racine Newton converge-t-il ? » est donc infiniment mal conditionne sur un ensemble de mesure nulle mais partout dense dans la frontiere.