import numpy as np
# Метод Ньютона для системы действительных уравнений.
# Комплексную задачу предварительно разделяют на действительные и мнимые составляющие.
def real_newton(residual, start, step=1e-4, tolerance=1e-7, maximum=20):
    # Создаём собственный массив неизвестных; начальное приближение сохраняется отдельно.
    x = np.array(start)
    for iteration in range(maximum):
        r = np.array(residual(x))
        # Критерий сходимости — максимальный модуль составляющей невязки.
        norm = max(abs(float(v)) for v in r)
        # Порядок вывода: номер итерации; максимальная невязка.
        print(iteration, norm)
        if norm <= tolerance:
            return x
        # Каждый столбец матрицы Якоби вычисляем центральной разностью.
        # Для каждой неизвестной возмущаем только одну координату.
        columns = []
        for k in range(len(x)):
            plus = np.array(x)
            minus = np.array(x)
            plus[k] += step
            minus[k] -= step
            rp = np.array(residual(plus))
            rm = np.array(residual(minus))
            columns.append((rp-rm)/(2*step))
        # Собираем матрицу по строкам уравнений и столбцам неизвестных.
        jacobian = np.array([[columns[k][j] for k in range(len(x))]
                             for j in range(len(x))])
        # Решаем систему J·Δx = −r для поправки к текущему состоянию.
        change = np.linalg.solve(jacobian, -r)
        # Уменьшаем длину шага, пока новая невязка не станет меньше прежней.
        # Пробное состояние принимаем только после независимого вычисления невязки.
        damping = 1.0
        accepted = False
        for trial in range(12):
            candidate = x + damping*change
            rc = residual(candidate)
            nc = max(abs(float(v)) for v in rc)
            if nc < norm:
                x = candidate
                accepted = True
                break
            damping *= 0.5
        # Отсутствие уменьшения невязки считается отказом, а не успешным завершением.
        if not accepted:
            raise ValueError("Демпфирование не уменьшило невязку")
    raise ValueError("Число итераций исчерпано")
# Аналитический эталон: два вещественных уравнения вместо комплексной производной.
# Проверочный пример содержит два независимых кубических уравнения.
# Здесь нет обращения к полевому решателю FEMM.
def control(x):
    return [x[0]+0.2*x[0]**3-2.0, x[1]+0.1*x[1]**3-1.0]
print(real_newton(control, [1.0, 0.5]))