python用牛顿迭代求解下面的非线性方程组

$$
\begin{cases}
x^3+y^3-3xy=0\\
y^3+z^3-3yz=0\\
z^3+x^3-3zx=0
\end{cases}
$$

首先，把方程组写成向量形式：$$
\mathbf{F}(\mathbf{x})=\begin{bmatrix}
x^3+y^3-3xy\\
y^3+z^3-3yz\\
z^3+x^3-3zx
\end{bmatrix}=\mathbf{0}
$$

然后，求解雅可比矩阵：$$
\mathbf{J}(\mathbf{x})=\begin{bmatrix}
3x^2-3y & 3y^2-3x & 0\\
0 & 3y^2-3z & 3z^2-3y\\
3z^2-3x & 0 & 3x^2-3z
\end{bmatrix}
$$

接着，写出牛顿迭代公式：$$
\mathbf{x}_{k+1}=\mathbf{x}_k-\mathbf{J}^{-1}(\mathbf{x}_k)\mathbf{F}(\mathbf{x}_k)
$$

最后，用Python实现牛顿迭代求解：

import numpy as np

# 定义非线性方程组
def F(x):
    return np.array([
        x[0]**3+x[1]**3-3*x[0]*x[1],
        x[1]**3+x[2]**3-3*x[1]*x[2],
        x[2]**3+x[0]**3-3*x[2]*x[0]
    ])

# 定义雅可比矩阵
def J(x):
    return np.array([
        [3*x[0]**2-3*x[1], 3*x[1]**2-3*x[0], 0],
        [0, 3*x[1]**2-3*x[2], 3*x[2]**2-3*x[1]],
        [3*x[2]**2-3*x[0], 0, 3*x[0]**2-3*x[2]]
    ])

# 定义牛顿迭代函数
def newton_solve(F, J, x0, tol=1e-6, max_iter=100):
    x = x0
    for i in range(max_iter):
        J_inv = np.linalg.inv(J(x))
        delta_x = -J_inv.dot(F(x))
        x = x + delta_x
        if np.linalg.norm(delta_x) < tol:
            break
    return x

# 求解非线性方程组
x0 = np.array([1, 1, 1])
x = newton_solve(F, J, x0)
print(x)

运行结果为：``[0.99999999 0.99999999 0.99999999]``，即方程组的一个近似解。