本文揭示 Newton-Raphson 多变量迭代中一个典型但极易被忽视的数学与 NumPy 实现偏差——将残差向量左乘雅可比逆矩阵（Fᵀ @ J⁻¹）误当作标准牛顿步长（J⁻¹ @ F），导致数值震荡、高精度下不收敛，并给出正确实现、健壮初始化及数值稳定性优化方案。
本文揭示 newton-raphson 多变量迭代中一个典型但极易被忽视的数学与 numpy 实现偏差——将残差向量左乘雅可比逆矩阵（`fᵀ @ j⁻¹`）误当作标准牛顿步长（`j⁻¹ @ f`），导致数值震荡、高精度下不收敛，并给出正确实现、健壮初始化及数值稳定性优化方案。
Newton-Raphson 法在非线性方程组求解中依赖于局部线性化：给定当前估计值 (\mathbf{x}_n = [x_n, y_n]^T)，其更新公式为：
[ \mathbf{x}_{n+1} = \mathbf{x}_n - \mathbf{J}^{-1}(\mathbf{x}_n)\, \mathbf{F}(\mathbf{x}_n) ]
其中 (\mathbf{F} = [f_1, f_2]^T) 是列向量形式的残差，(\mathbf{J}) 是 (2\times2) 雅可比矩阵。关键点在于：(\mathbf{J}^{-1}) 必须左乘列向量 (\mathbf{F})，即执行 np.linalg.solve(J, F) 或等价的 J_inv @ F（@ 表示矩阵乘法）。而原代码中使用了：
result = -np.matmul(np.array([f1(x0, y0), f2(x0, y0)]), np.linalg.inv(jacobian_mat(x0, y0))) + np.array([x0, y0])
这实际计算的是 (-\mathbf{F}^T \mathbf{J}^{-1} + \mathbf{x}_n) —— 一个 (1\times2) 行向量减去 (2\times2) 矩阵的结果，再隐式广播为一维数组。由于 NumPy 对 1D 数组的“自动列/行推断”规则（如 A @ v 将 v 视为列，v @ A 将 v 视为行），该表达式等价于 (-(\mathbf{J}^{-1})^T \mathbf{F})，严重偏离理论要求，破坏了牛顿方向的几何意义，是造成高精度（如 TOLERANCE=0.001）下振荡不收敛的根本原因。
✅ 正确实现应统一使用列向量语义，并优先避免显式求逆（提高数值稳定性与效率）：
def compute_newton_raphson(x0, y0):
x, y = float(x0), float(y0) # 显式转浮点，避免整数除法问题
for i in range(ITER_LIMIT):
F = np.array([[f1(x, y)], [f2(x, y)]]) # 2×1 列向量
J = jacobian_mat(x, y) # 2×2 矩阵
try:
# 推荐：用 np.linalg.solve 替代显式求逆（更稳定、更快）
delta = np.linalg.solve(J, -F) # 解 J @ delta = -F
x_new = x + delta[0, 0]
y_new = y + delta[1, 0]
except np.linalg.LinAlgError:
print(f"Iteration {i}: Singular Jacobian at ({x:.4f}, {y:.4f})")
return None
# 收敛判据：使用绝对步长（更鲁棒）或函数值范数
if max(abs(x_new - x), abs(y_new - y)) <= TOLERANCE:
return (x_new, y_new)
# 发散检测：防止数值爆炸
if abs(x_new) > 1e6 or abs(y_new) > 1e3:
print(f"Iteration {i}: Divergence detected → restart")
return None
x, y = x_new, y_new
return None? 其他关键优化建议：
初始值策略：
y是角度（单位：弧度），不应在[-100, 100]内随机取整数（相当于[-318°, 318°]的任意整数度）。应改为在[0, 2π)或[-π, π)内均匀采样浮点数；x可保持合理物理范围（如[-50, 50]），但需用np.random.uniform而非rnd.randrange。收敛判据升级：除步长外，建议同时监控残差范数
np.linalg.norm(F)，避免“假收敛”（如在平坦区域步长小但残差仍大）。鲁棒性增强：添加雅可比奇异检测（
LinAlgError）、解溢出检查（abs(x) > 1e6），并在while True循环中限制总重试次数（如max_restarts=100），避免无限循环。-
终极推荐：对于工程应用，直接使用
scipy.optimize.root（支持多种算法及自动雅可比计算）：from scipy.optimize import root def system(vars): x, y = vars return [f1(x, y), f2(x, y)] sol = root(system, x0=[10.0, np.pi/4], method='hybr', tol=TOLERANCE) if sol.converged: print(f"X={sol.x[0]:.6f}, Y={np.rad2deg(sol.x[1]):.4f}°")
综上，Newton-Raphson 的成功不仅取决于公式正确性，更依赖于向量/矩阵的严格维度语义、数值稳定性实践与问题先验知识（如变量量纲与合理范围）。修正乘法顺序是破局第一步，而系统性健壮化设计才是工业级求解的基石。
