求解非线性方程组的Newton迭代法和逆Broyden迭代法+python

import numpy as np


# 定义非线性方程组
def f(x):
    f1 = x[0] - 0.7 * np.sin(x[0]) - 0.2 * np.cos(x[1])
    f2 = x[1] - 0.7 * np.cos(x[0]) + 0.2 * np.sin(x[1])
    return np.array([f1, f2])


# 定义雅可比矩阵
def jacobian(x):
    df1_dx1 = 1 - 0.7 * np.cos(x[0])
    df1_dx2 = 0.2 * np.sin(x[1])
    df2_dx1 = 0.7 * np.sin(x[0])
    df2_dx2 = 1 + 0.2 * np.cos(x[1])
    return np.array([[df1_dx1, df1_dx2], [df2_dx1, df2_dx2]])


# Newton 迭代法
def newton_method(x0, max_iter=3):
    x = x0.copy()
    for k in range(max_iter):
        print(f"\n=== Newton 第 {k + 1} 次迭代 ===")
        print("当前解 x =", x)
        F = f(x)
        print("F(x) =", F)
        J = jacobian(x)
        print("雅可比矩阵 J =\n", J)

        # 解线性方程组 J * dx = -F
        delta_x = np.linalg.solve(J, -F)
        print("步长 Δx =", delta_x)

        # 更新解
        x += delta_x
        print("更新后 x =", x)

    return x


# 逆 Broyden 迭代法
def inverse_broyden_method(x0, max_iter=3):
    x = x0.copy()
    n = len(x)
    B = np.eye(n)  # 初始 B 矩阵为单位矩阵

    for k in range(max_iter):
        print(f"\n=== Broyden 第 {k + 1} 次迭代 ===")
        print("当前解 x =", x)
        F = f(x)
        print("F(x) =", F)

        # 计算搜索方向 d = -B * F
        d = -B @ F
        print("搜索方向 d =", d)

        # 保存旧值
        x_prev = x.copy()
        F_prev = F.copy()

        # 更新解
        x += d
        print("更新后 x =", x)

        # 计算 s 和 y
        s = x - x_prev
        F_new = f(x)
        y = F_new - F_prev
        print("s =", s)
        print("y =", y)

        # Broyden 更新公式
        # numerator = np.outer((y - B @ s), s.T @ B)
        numerator = np.dot((y - B @ s), s.T @ B)
        denominator = s.T @ B @ s
        B += numerator / denominator
        print("更新后的 B 矩阵 =\n", B)

    return x


# 主程序
if __name__ == "__main__":
    x0 = np.array([0.5, 0.5])
    print("\n=== Newton 迭代法结果 ===")
    newton_solution = newton_method(x0)
    print("\n最终解 (Newton):", newton_solution)
    print("F(x) =", f(newton_solution))

    print("\n=== Broyden 迭代法结果 ===")
    broyden_solution = inverse_broyden_method(x0)
    print("\n最终解 (Broyden):", broyden_solution)
    print("F(x) =", f(broyden_solution))