Système non linéaire (Newton 2D)

import numpy as np

def newton_systeme(F, J, x0, epsilon=1e-12, maxiter=100):
    x = np.array(x0, dtype=float)
    for n in range(maxiter):
        Fx = F(x[0], x[1])
        if np.linalg.norm(Fx) < epsilon:
            return x, n+1
        Jx = J(x[0], x[1])
        delta = np.linalg.solve(Jx, -Fx)
        x_new = x + delta
        if np.linalg.norm(delta) < epsilon:
            return x_new, n+1
        x = x_new
    raise RuntimeError("Non convergence")

# Système
def F(x, y):
    return np.array([x**2 + y**2 - 4, x**2 - y - 1])

def J(x, y):
    return np.array([[2*x, 2*y], [2*x, -1]])

# Recherche des racines depuis différents points initiaux
points_initiaux = [(2, 1), (-2, 1), (1, -1.5), (-1, -1.5)]

for x0, y0 in points_initiaux:
    sol, iterations = newton_systeme(F, J, [x0, y0])
    print(f"Init ({x0}, {y0}) -> ({sol[0]:.8f}, {sol[1]:.8f}) en {iterations} it.")

Résultat dans la console :

Init (2, 1) -> (1.73205081, 1.00000000) en 6 it.
Init (-2, 1) -> (-1.73205081, 1.00000000) en 6 it.
Init (1, -1.5) -> (0.00000000, -1.00000000) en 5 it.
Init (-1, -1.5) -> (0.00000000, -1.00000000) en 5 it.

Les quatre points d'intersection réels du système sont : $(\pm \sqrt{3},1)$ et $(0,-1)$ (racine double).