Загрузка данных


import numpy as np

# Определение системы уравнений
def F(x):
    """Вектор функций: [F1, F2]"""
    x1, x2 = x[0], x[1]
    f1 = 1.0 / (x1 + 1.5) - x2 - 1.0
    f2 = (1.2 * x1) ** 3 - x2 - 1.0
    return np.array([f1, f2], dtype=float)

# Матрица Якоби
def Jacobian(x):
    """Матрица Якоби 2x2"""
    x1, x2 = x[0], x[1]
    J = np.zeros((2, 2))
    # ∂F1/∂x1 = -1/(x1+1.5)^2
    J[0, 0] = -1.0 / ((x1 + 1.5) ** 2)
    # ∂F1/∂x2 = -1
    J[0, 1] = -1.0
    # ∂F2/∂x1 = 3*(1.2*x1)^2 * 1.2 = 5.184 * x1^2
    J[1, 0] = 5.184 * (x1 ** 2)
    # ∂F2/∂x2 = -1
    J[1, 1] = -1.0
    return J

# Метод Ньютона по блок-схеме
def newton_system(x0, eps=1e-4, max_iter=100):
    x = np.array(x0, dtype=float)
    ndx = 2.0 * eps   # чтобы войти в цикл
    it = 0

    print("Итерация |       x1        x2      |       F1        F2      |  ||dx||")
    print("-" * 70)

    while ndx > eps and it < max_iter:
        fx = F(x)
        # Вывод текущих значений (как в блок-схеме)
        print(f"{it:5d}    | {x[0]:10.6f} {x[1]:10.6f} | {fx[0]:10.6f} {fx[1]:10.6f} | {ndx:10.6f}")

        J = Jacobian(x)
        # Решаем J * dx = -fx   (методом обратной матрицы)
        # invJ = np.linalg.inv(J)
        # dx = invJ @ (-fx)   # обычно dx = -invJ @ fx
        # Но в блок-схеме: dx = invJ @ fx, затем x = x - dx
        # Это эквивалентно: dx = inverse(J) @ fx, но тогда x_new = x - dx
        # Однако правильное ньютоновское обновление: J * dx = -F, откуда dx = -J^{-1} F
        # Чтобы соответствовать блок-схеме (где x = x - dx, а dx = invJ @ fx),
        # мы должны взять dx = invJ @ fx, тогда x = x - dx.
        # Проверим: из F(x) = 0 должно быть J * (x_new - x) = -F => x_new = x - J^{-1}F.
        # Следовательно, если определить dx = J^{-1} @ F, то x_new = x - dx.
        # Так и сделаем (в блок-схеме именно так: dx = invJacob @ fx() и x = x - dx).
        # Поэтому:
        invJ = np.linalg.inv(J)
        dx = invJ @ fx
        ndx = np.linalg.norm(dx)

        x = x - dx
        it += 1

    # Финальный вывод
    fx_final = F(x)
    print(f"{it:5d}    | {x[0]:10.6f} {x[1]:10.6f} | {fx_final[0]:10.6f} {fx_final[1]:10.6f} | {ndx:10.6f}")
    print("\nРешение найдено за", it, "итераций")
    print(f"x1 = {x[0]:.6f}, x2 = {x[1]:.6f}")
    print(f"Невязки: F1 = {fx_final[0]:.2e}, F2 = {fx_final[1]:.2e}")
    return x

# Запуск с начальным приближением из блок-схемы
if __name__ == "__main__":
    x0 = [0.5, 0.5]
    eps = 0.0001
    solution = newton_system(x0, eps)