import numpy as np

# Те же F(x) из предыдущего кода
def F(X):
    x, y, z = X
    return np.array([
        x**2 + x - y*z - 0.1,
        -3*y + y**2 - x*z + 0.2,
        2*z + z**2 + x*y - 0.3
    ], dtype=float)

# Якобиан J(x)
def J(X):
    x, y, z = X
    return np.array([
        [2*x + 1,   -z,       -y],
        [  -z,    2*y - 3,    -x],
        [   y,       x,    2 + 2*z]
    ], dtype=float)

# Твой Гаусс
def gauss_solve(A, bb):
    n = A.shape[0]
    M = np.hstack((A.copy(), bb.reshape(-1, 1)))

    # Прямой ход с выбором главного элемента по столбцу
    for i in range(n):
        max_row = i
        for k in range(i + 1, n):
            if abs(M[k, i]) > abs(M[max_row, i]):
                max_row = k
        M[[i, max_row]] = M[[max_row, i]]  # перестановка строк

        for k in range(i + 1, n):
            if M[i, i] == 0:
                raise ValueError("Нулевой ведущий элемент в Гауссе")
            c = M[k, i] / M[i, i]
            M[k, i:] -= c * M[i, i:]

    # Обратный ход
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = M[i, n]
        for j in range(i + 1, n):
            x[i] -= M[i, j] * x[j]
        x[i] /= M[i, i]

    return x

def gauss_inverse(A):
    n = A.shape[0]
    A_inv = np.zeros((n, n))
    I = np.eye(n)

    for i in range(n):
        e = I[:, i]
        A_inv[:, i] = gauss_solve(A, e)

    return A_inv

def newton(x0, eps=1e-6, max_iter=10):
    X = np.array(x0, dtype=float)

    print("\nМетод Ньютона (ε = %.0e)" % eps)
    print("k\t     x_k\t\t     y_k\t\t     z_k\t\tmax|F_i(x_k)|")

    for k in range(1, max_iter + 1):
        Fx = F(X)
        Jx = J(X)
        Jinv = gauss_inverse(Jx)      # обратная к Якоби
        delta = -Jinv.dot(Fx)
        X_new = X + delta
        r_new = F(X_new)
        r_max = np.max(np.abs(r_new))

        print(f"{k:2d}\t{X_new[0]: .8f}\t{X_new[1]: .8f}\t"
              f"{X_new[2]: .8f}\t{r_max: .2e}")
        print("J(x_k) =\n", Jx)
        print("J(x_k)^(-1) =\n", Jinv, "\n")

        # критерий остановки по шагу
        if np.max(np.abs(delta)) < eps:
            print("Достигнута точность по шагу ε =", eps)
            break

        X = X_new

    return X_new

# Тот же старт
x0 = [0.10, 0.07, 0.15]
root_newton = newton(x0)

print("\nИтог (метод Ньютона):")
print("x =", root_newton[0])
print("y =", root_newton[1])
print("z =", root_newton[2])
print("Невязки F(x):", F(root_newton))