迭代求解非线性方程组的代码

下面是一个简单的迭代求解非线性方程组的代码，使用的是牛顿迭代法：

import numpy as np

def newton(func, x0, tol=1e-6, max_iter=100):
    """
    用牛顿迭代法求解非线性方程组的解

    :param func: 函数f(x)，输入x是一个向量，输出也是一个向量
    :param x0: 初始解，可以是一个向量或者一个标量
    :param tol: 精度要求
    :param max_iter: 最大迭代次数
    :return: 返回解向量
    """
    x = x0
    for i in range(max_iter):
        fx = func(x)
        J = jacobian(func, x)
        dx = np.linalg.solve(J, -fx)
        x = x + dx
        if np.linalg.norm(dx) < tol:
            break
    return x

def jacobian(func, x, eps=1e-6):
    """
    计算函数f(x)在x处的雅可比矩阵

    :param func: 函数f(x)，输入x是一个向量，输出也是一个向量
    :param x: 输入向量
    :param eps: 数值微分的精度
    :return: 雅可比矩阵
    """
    n = len(x)
    J = np.zeros((n, n))
    for i in range(n):
        x1 = np.array(x)
        x1[i] += eps
        fx1 = func(x1)
        x2 = np.array(x)
        x2[i] -= eps
        fx2 = func(x2)
        J[:, i] = (fx1 - fx2) / (2 * eps)
    return J

使用示例：

假设我们要求解方程组：

$$\begin{cases} x^2+y^2=1 \\ x^3-y=0 \end{cases}$$

定义函数：

def func(x):
    return np.array([x[0]**2 + x[1]**2 - 1, x[0]**3 - x[1]])

调用 `newton` 函数：

x0 = [1, 1]
x = newton(func, x0)
print("解向量：", x)

输出结果：

解向量： [0.71414284 0.70014004]