本文深入剖析 newton-raphson 法在求解二维向量环方程组时收敛失效的根本原因——核心在于雅可比矩阵逆运算与残差向量乘法顺序错误，并提供符合数学定义的正确实现、鲁棒初始化策略及数值稳定性优化建议。
本文深入剖析 newton-raphson 法在求解二维向量环方程组时收敛失效的根本原因——核心在于雅可比矩阵逆运算与残差向量乘法顺序错误，并提供符合数学定义的正确实现、鲁棒初始化策略及数值稳定性优化建议。
Newton-Raphson 方法求解非线性方程组 $ \mathbf{F}(\mathbf{x}) = \mathbf{0} $ 的标准迭代格式为：
$$ \mathbf{x}_{n+1} = \mathbf{x}_n - \mathbf{J}^{-1}(\mathbf{x}_n)\, \mathbf{F}(\mathbf{x}_n) $$
其中 $ \mathbf{J} $ 是雅可比矩阵（2×2），$ \mathbf{F} = [f_1, f_2]^\top $ 是列向量形式的残差。关键点在于：$ \mathbf{J}^{-1} $ 必须左乘（作用于）列向量 $ \mathbf{F} $，即 np.linalg.solve(J, F) 或等价的 np.linalg.inv(J) @ F。
而原代码中使用了：
result = -np.matmul(np.array([f1(x0, y0), f2(x0, y0)]), np.linalg.inv(jacobian_mat(x0, y0))) + np.array([x0, y0])
这等价于 $ \mathbf{x}_{n+1} = \mathbf{x}_n - \mathbf{F}^\top \mathbf{J}^{-1} $，即用行向量右乘逆雅可比矩阵——数学上无定义（维度不匹配），且实际计算中因 NumPy 对一维数组的隐式广播规则，等效于 $ \mathbf{x}_n - (\mathbf{J}^{-\top}\mathbf{F})^\top $，引入了转置误差。该错误导致迭代方向失真，即使在粗精度（如 TOLERANCE=0.01）下偶然“凑效”，在更高精度要求（如 0.001）时必然暴露发散或振荡行为。
✅ 正确实现应严格遵循列向量约定。以下是修复后的核心函数：
def compute_newton_raphson(x0, y0):
x, y = x0, y0
for i in range(ITER_LIMIT):
F = np.array([f1(x, y), f2(x, y)]) # 列向量 [f1; f2]
J = jacobian_mat(x, y) # 2x2 矩阵
# ✅ 推荐：用 solve 替代 inv（更稳定、高效）
try:
delta = np.linalg.solve(J, F) # 解 J @ delta = F → delta = J^{-1} F
except np.linalg.LinAlgError:
print(f"Iteration {i}: Singular Jacobian at ({x:.4f}, {y:.4f})")
return None
x_new = x - delta[0] # x_{n+1} = x_n - (J^{-1}F)[0]
y_new = y - delta[1] # y_{n+1} = y_n - (J^{-1}F)[1]
# 使用残差范数或步长范数判断收敛（更鲁棒）
if np.linalg.norm([x_new - x, y_new - y], ord=np.inf) <= TOLERANCE:
return (x_new, y_new)
# 防发散：检测变量是否失控增长
if abs(x_new) > 1e6 or abs(y_new) > 1e3:
print(f"Iteration {i}: Divergence detected at ({x_new:.4f}, {y_new:.4f})")
return None
x, y = x_new, y_new
print(f"Failed to converge within {ITER_LIMIT} iterations.")
return None⚠️ 其他关键改进建议：
-
初始化策略优化：
原代码对角度y在[-100, 100]弧度范围内随机取整，极不合理（100 rad ≈ 16π，覆盖约 16 个完整周期）。应改为在物理合理区间（如[-π, π]或[-2π, 2π]）内均匀采样浮点数：y0 = np.random.uniform(-np.pi, np.pi) # 角度初始化更合理 x0 = np.random.uniform(-50, 50) # 位移量按工程量纲设定
-
收敛判据升级：
仅检查步长|x_{n+1}-x_n|和|y_{n+1}-y_n|不足以保证解的精度。强烈建议同时监控残差范数：if np.linalg.norm(F, ord=np.inf) < TOLERANCE * 1e-2: # 残差足够小 return (x, y) 健壮性增强：
在while True循环中增加最大重试次数（如max_retries = 1000），避免无限等待；对奇异雅可比矩阵（LinAlgError）和数值溢出做显式捕获与跳过。
总结：Newton-Raphson 的成功高度依赖于数学公式严格实现与数值实践合理性。修正矩阵乘法顺序是解决收敛失效的前提；辅以物理意义明确的初值、多重收敛判据及异常处理机制，才能构建稳定、高精度的非线性求解器。对于工程问题，亦可直接采用 scipy.optimize.root(method='hybr') 等成熟封装，降低实现风险。
